Paper deep dive
A Hierarchical Multi-Agent System for Autonomous Discovery in Geoscientific Data Archives
Dmitrii Pantiukhin, Ivan Kuznetsov, Boris Shapkin, Antonia Anna Jost, Thomas Jung, Nikolay Koldunov
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 7/20/2026, 1:07:17 PM
Summary
The paper introduces PANGAEA-GPT, a hierarchical multi-agent system designed for autonomous data discovery and analysis in geoscientific data archives like PANGAEA. The system utilizes a Supervisor-Worker topology where a Supervisor Agent decomposes natural language queries into subtasks and routes them to specialist worker agents (Oceanographer, Ecologist, Visualization, DataFrame, Writer) for execution in a sandboxed environment. The framework emphasizes data-type-aware routing, self-correction via execution feedback, and reflexive quality control. The authors validate the system through four scenarios involving microplastic distribution, reanalysis validation, trajectory analysis, and biodiversity coupling, demonstrating its ability to handle complex, multi-step workflows with minimal human intervention. Additionally, a search benchmark shows that the agentic search approach significantly outperforms baseline keyword matching and simple LLM translation in precision and parameter coverage.
Entities (13)
Relation Signals (15)
Agentic Search → outperforms → Simple LLM
confidence 95% · compared to the Simple LLM’s 5.46/10
Agentic Search → outperforms → Elasticsearch BM25
confidence 95% · The Agentic Search strategy increased the overall mean score to 8.14/10 compared to ... the Baseline’s 2.87/10
PANGAEA-GPT → usestopology → Supervisor-Worker topology
confidence 95% · our architecture implements a centralized Supervisor-Worker topology
Supervisor Agent → delegatesto → DataFrame Agent
confidence 90% · routes each subtask to one of five specialist worker agents: ... (iv) a DataFrame Agent
Supervisor Agent → delegatesto → Writer Agent
confidence 90% · routes each subtask to one of five specialist worker agents: ... (v) a Writer Agent
Supervisor Agent → delegatesto → Oceanographer Agent
confidence 90% · The Supervisor then routes each subtask to one of five specialist worker agents: (i) an Oceanographer Agent
Supervisor Agent → delegatesto → Ecologist Agent
confidence 90% · routes each subtask to one of five specialist worker agents: ... (ii) an Ecologist Agent
Supervisor Agent → delegatesto → Visualization Agent
confidence 90% · routes each subtask to one of five specialist worker agents: ... (iii) a Visualization Agent
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The rapid accumulation of Earth science data has created a significant scalability challenge; while repositories like PANGAEA host vast collections of datasets, citation metrics indicate that a substantial portion remains underutilized, limiting data reusability. Here we present PANGAEA-GPT, a hierarchical multi-agent framework designed for autonomous data discovery and analysis. Unlike standard Large Language Model (LLM) wrappers, our architecture implements a centralized Supervisor-Worker topology with strict data-type-aware routing, sandboxed deterministic code execution, and self-correction via execution feedback, enabling agents to diagnose and resolve runtime errors. Through use-case scenarios spanning physical oceanography and ecology, we demonstrate the system's capacity to execute complex, multi-step workflows with minimal human intervention. This framework provides a methodology for querying and analyzing heterogeneous repository data through coordinated agent workflows.
Tags
Links
- Source: https://arxiv.org/abs/2602.21351v1
- Canonical: https://arxiv.org/abs/2602.21351v1
Trouble viewing inline? Open PDF directly →
Full Text
95,080 characters extracted from source content.
Expand or collapse full text
A Hierarchical Multi-Agent System for Autonomous Discovery in Geoscienti c Data Archives Dmitrii Pantiukhin 1 , Ivan Kuznetsov 1 , Boris Shapkin 1 , Antonia Anna Jost 1 , Thomas Jung 1 , and Nikolay Koldunov 1 1 Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany Abstract The rapid accumulation of Earth science data has created a signi cant scalability challenge; while repositories like PANGAEA host vast collections of datasets, citation metrics indicate that a substantial portion remains underutilized, limiting data reusabil- ity. Here we present PANGAEA-GPT, a hierarchical multi-agent frameworkdesigned for autonomous data discovery and analysis. Unlike standard Large Language Model (LLM) wrappers, our architecture implements a centralized Supervisor-Worker topology with strict data-type-aware routing, sandboxed deterministic code execution, and self- correction via execution feedback, enabling agents to diagnose and resolve runtime errors. Through use-case scenarios spanning physical oceanography and ecology, we demonstrate the system’s capacity to execute complex, multi-step work ows with minimal human intervention. This framework provides a methodology for querying and analyzing hetero- geneous repository data through coordinated agent work ows. 1 Introduction The increasing volume of digital data from sources such as satellite observations, autonomous underwater vehicles, and research vessels has increased the demand for scalable, automated analysis methods. Major international repositories, including PANGAEA, NASA’s Earth Ob- serving System Data and Information System (EOSDIS), and the Copernicus Climate Data Store (C3S), now support essential climate and oceanographic research, collectively hosting petabytes of observational and model-derived data [ Kobler et al.,1995,Felden et al.,2023, Li et al.,2023b]. PANGAEA alone catalogs over 400,000 curated datasets spanning more than 800 distinct geoscienti c parameters, from Conductivity-Temperature-Depth (CTD) hy- drographic pro les and biogeochemical uxes to ice-core proxy records and coupled model ensembles [ Felden et al.,2023]. Yet the scale of curation exceeds the rate of scienti c reuse: analysis of the Data Citation Index reveals that nearly 90% of archived datasets have never been cited in the peer-reviewed literature [ Robinson-Garc’ıa et al.,2016]. This usability gap stems from interoperability challengesacross the data lifecycle. In- complete or inconsistent documentation obscures essential context (measurement provenance, instrument calibration history, quality-control ags), without which downstream reuse is im- practical [ Gil et al.,2016]. Format heterogeneity ampli es the problem: the same physical quantity may be encoded as a column in a at CSV, a variable in a multi-dimensional NetCDF array, or a chunked Zarr store, each demanding a distinct software toolchain and non-trivial 1 arXiv:2602.21351v1 [cs.AI] 24 Feb 2026 domain expertise to interpret correctly [Gil et al.,2016,Schnase et al.,2016]. While the FAIR (Findable, Accessible, Interoperable, Reusable) data principles provide an essential framework for addressing these barriers [ Wilkinson et al.,2016], implementing them at scale remains chal- lenging even for dedicated community e orts [ Marsden et al.,2025]. The cumulative technical load required to discover, harmonize, and analyze disparate datasets remains a primary bottle- neck to data reuse [ Li et al.,2023b]. Large Language Models (LLMs) have evolved from text generators into autonomous agents that reason over problems, decompose tasks, and invoke external tools [ Boiko et al.,2023a, Schick et al.,2023]. This shift has given rise to Multi-Agent Systems (MAS) in which com- plex problems are partitioned across specialized agents [ Hong et al.,2024,Li et al.,2023a, Guo et al.,2024b]. Hierarchical agent architectures have demonstrated the ability to navi- gate knowledge graphs and design experiments in materials science [ Ghafarollahi and Buehler, 2024b,a]. Similarly, the AI Scientist-v2 [Yamada et al.,2025] and Google’s AI Co-Scientist [ Gottweis et al.,2025] have been applied to automating research cycles and hypothesis gener- ation, while AlphaEvolve [ Novikov et al.,2025] uses LLM-driven evolutionary computation to discover and optimize algorithms, advancing results in mathematics and computing infrastruc- ture. In laboratory settings, systems like Coscientist [ Boiko et al.,2023b] and ChemCrow [Bran et al. ,2024] have achieved autonomous planning in chemical synthesis.More recently, Robin [ Ghareeb et al.,2025] has automated hypothesis-driven biomedical discovery, Data Interpreter [ Hong et al.,2025] targets general-purpose data science tasks, and DS-Agent [Guo et al.,2024a] automates machine learning pipelines. Current LLMs already rival human domain experts on graduate-level science benchmarks [ Rein et al.,2024], and their expanding context windows enable coherent multi-step tool use. Orchestration frameworks such as HuggingGPT show that an LLMcontroller can dynamically route subtasks to specialist models [ Shen et al.,2023]. Whether these capabilities transfer to the messier reality of scienti c data archives, with their inconsistent metadata and format heterogeneity, remains an open question. Within geosciences, LLM adoption has progressed from domain-speci c foundation mod- els to tool-augmented autonomous agents. Models like GeoGalactica [ Lin et al.,2023] and OceanGPT [ Bi et al.,2024] established domain-speci c linguistic capabilities. Inthe domain of climate services, ClimSight [ Koldunov and Jung,2024,Kuznetsov et al.,2025] demonstrated that augmenting LLMs with localized climate model output, RAG-based report retrieval, and an agent-based architecture can deliver actionable, location-speci c climate assessments at scale. Tool-augmented models like GeoGPT [ Zhang et al.,2023b] enabled geospatial data processing, followed by single-agent frameworks like GeoAgent [ Chen et al.,2024] and LLM-Find [Ning et al. ,2025]. However, as identi ed in our prior architectural survey [Pantiukhin et al.,2025], none of these systems provide end-to-end integration with the retrieval, format-handling, and analytical infrastructure of a production-scale geoscienti c data repository. Despite this progress, existing approaches address individual components of the scienti c work ow in isolation (discovery, analysis, or visualization) without integrating them into a uni- ed, self-correcting pipeline capable of operating acrossthe format diversity of a major reposi- tory. Moreover, no multi-agent system has been deployed forgeoscienti c data archives, where the inherent disciplinary breadth, spanning physical oceanography, biogeochemistry, paleocli- matology, and ecology, demands specialist agents with distinct analytical toolchains. What is missing is an architecture that combines semantic search, data-type-aware task routing, and deterministic code execution in a single orchestrated framework. Here, we describe PANGAEA-GPT, a hierarchical multi-agent system rst outlined in Pan- tiukhin et al.[2025], and provide its full architecture along with the rst scenario-driven evalu- ation on real research work ows. We demonstrate that by combining an iterative, self-re ning agentic search with a supervisor-orchestrated team of specialist data agents, the system can 2 execute complex, multi-step Earth science analyses with reduced manual intervention. 2 Results 2.1 System overview The system operates through a hierarchical Supervisor-Worker topology that enforces a strict separation between reasoning and execution (Fig.1). A natural language query is received directly by the Supervisor Agent, which decomposes it into subtasks based on the schema, dimensionality, and coordinate reference system of each retrieved dataset. The Supervisor then routes each subtask to one of ve specialist worker agents: (i) an Oceanographer Agent for N-dimensional arrays (NetCDF, Zarr) and external climate reanalysis APIs; (i) an Ecologist Agent for biodiversity metrics and taxonomic aggregation; (i) a Visualization Agent with a re exive quality-control loop and a retrieval-augmentedgeneration (RAG) index of domain- speci c exemplar plots; (iv) a DataFrame Agent for statistical computation on tabular data; and (v) a Writer Agent for narrative synthesis via a vision-language model interface (Methods 5.2.2). This modality-based routing ensures that, for example, gridded arrays are never processed by tabular agents. 2.2 Validation scenarios We evaluated PANGAEA-GPT on four scenarios that collectively test the full scienti c work- ow: data retrieval, cross-domain integration, statistical analysis, and visualization. 2.2.1 Scenario 1: Microplastic distribution in the Weddell Sea Quantifying the coupling between ocean circulation and plastic accumulation typically requires bridging two distinct data regimes: Lagrangian point observations and Eulerian ow elds. To test the system’s capacity for this cross-domain integration, we loaded a published dataset of sea-surface microplastic concentrations from the WeddellSea (N= 34 Manta trawl stations; DOI: 10.1594/PANGAEA.941076) and directed the agents to integrate it with ocean reanalysis elds. Guided by the user’s de nition of the study period, the Oceanographer Agent resolved the speci c spatiotemporal constraints, retrieved daily surface velocity elds from the Coperni- cus Marine Service (GLORYS12V1), and computed a temporal meancirculation eld over the months of the expedition. The Visualization Agent rendered a regional map (Fig. 2a) overlaying current speed, ow streamlines, and stations scaled by microplastic abundance. The generated streamlines visually captured the cyclonic rotation of theWeddell Gyre. During execution, the agents autonomously resolved one API-level error and one visualization layout failure (Sup- plementary Table S1). The system also handled domain-speci c constraints of the underlying ocean model. When the initial request for surface-level datawas rejected by the service, which de nes its shallowest vertical layer at 0.494 m rather than the sea surface, the Oceanographer Agent autonomously parsed the error response, identi ed thenearest valid vertical coordinate, and resubmitted the retrieval without human intervention.In the visualization work ow, the re exive critic identi ed that the legend obscured the central streamlines, prompting the agent to refactor the gure layout into dedicated colorbar and legend axes. Following a user request to quantify the coupling, the DataFrame Agent performed a sta- tistical evaluation onN= 21 co-located stations (positive abundance, valid current speed), revealing a moderate negative association between currentspeed and microplastic concentra- tion (Spearmanρ=−0.47,p= 0.033; Fig. 2b). A Mann–WhitneyUtest con rmed signi - 3 Phase 1: Data Discovery User Query (Natural Language) Search Agent Search Modes Simple: Optimized Boolean Query Deep: Parallel Queries & Consolidation PANGAEA Archive Publication Data (CrossRef) Selected Datasets & Analysis Request Phase 2: Data Analysis & Visualization Supervisor AgentPlanning Tool Stateful Workflow (LangGraph) Data-Type Aware Routing Task Delegation OceanographerEcologistVisualizationDataFrame Isolated Execution Environment (Sandbox) Toolboxes Core Tools (Shared) Python REPL (Jupyter K.)File Lister Code Examples (RAG)Package Installer Wise Agent (Advisor)Reflection (Vision) Validator (Vision Model) *DataFrame Agent uses only REPL & File Lister Oceanographer Specific Tools ERA5 RetrievalCopernicus Marine Synthesis Writer Agent Describe Images Tool Scientific Report & Visualizations Figure 1: Conceptual framework of the PANGAEA-GPT Multi-Agent System (MAS). The system uses a two-phase, hierarchical architecture. First,a Search Agent discovers relevant datasets from a user’s natural language request. A Supervisor Agent then delegates analysis and visualization tasks to a team of specialist agents working within a secure sandbox environment. Finally, a Writer Agent synthesizes the results into a cohesive report. cantly higher abundance in the low-speed regime (U= 88.0,p= 0.022; median 2,558 vs. 268 items km −2 ). We emphasize that these statistics serve to demonstrate the system’s capacity to autonomously execute a complete analytical pipeline, from data retrieval through statisti- cal testing, rather than to draw robust ecological conclusions. With only 21 stations and no control for confounding variables, these results are exploratory and should not be treated as con rmatory. 2.2.2 Scenario 2: Reanalysis validation against mooring observations Validating global climate models against in-situ observations typically requires extensive manual coding to bridge the gap between disparate tabular data and multi-dimensional gridded outputs. We tasked the system with validating temperature records from ten mooring datasets deployed at the LTER HAUSGARTEN observatory (Fram Strait, 2006–2012; PANGAEA DOIs 845608– 4 (a) Mean surface current speed (January–February 2018– 2019, GLORYS12V1) with streamlines and Manta trawl sta- tions scaled by microplastic abundance. The visualization was iteratively re ned via the re exive quality control loop (3/10→9/10). (b) Log-linear regression of abun- dance against current speed (N= 21). Annotations: Spearmanρ= −0.47 (p= 0.033); Mann–Whitney p= 0.022. Statistical annota- tions were automatically generated and placed on the plot by the Visu- alization Agent. Figure 2:Exploratory analysis of microplastic distribution in the Weddell Sea(Sup- plementary Note 1). 845617) against the Copernicus Marine Service global reanalysis (GLORYS12V1). The work ow proceeded through three stages. First, the Visualization Agent generated exploratory depth-coverage plots. Second, the DataFrame Agent programmatically harmo- nized the ten CSV les, standardizing column names and ltering for upper-ocean records (N= 135,678). Finally, the Oceanographer Agent performed the 4D matchup. During this step, it encountered a type mismatch between the mooring timestamps and the reanalysis model’s native date format. The agent diagnosed the incompatibility, refactored the temporal alignment, and re-executed the matchup, completing all 135,678 co-locations without human intervention (Supplementary Note 2). The resulting validation scatter plot (Fig. 3) shows a warm bias of +0.35 ◦ C, root-mean-square error (RMSE) of 1.09 ◦ C, and Pearsonr= 0.31. The positive bias is consistent with known o sets in the GLORYS12V1 reanalysis that are ampli ed in regions of intense Atlantic Water in ow such as the Fram Strait [ Lellouche et al.,2021], and the moderate correlation re ects the use of nearest-grid-point matching without spatial inter- polation across steep lateral temperature gradients between Atlantic and Arctic water masses. 2.2.3 Scenario 3: Resource-ecient retrieval and trajectory analysis Validating gridded reanalysis products against moving platforms presents challenges regard- ing data transfer eciency and geometric mismatch. We tasked the system with comparing ERA5 2-meter air temperature against shipboard measurements along the RV Polarstern track during the MOSAiC expedition (PS122/1, September–December2019; DOI: 10.1594/PAN- GAEA.935221 ). The DataFrame Agent parsed the navigation log to extract spatiotemporal bounds, which the Oceanographer Agent used to execute a targeted cloud-native retrieval from the Arraylake ERA5 data store [ Earthmover,2025]. The agent then interpolated the gridded eld to the ship’s 2,013 hourly positions, producing a co-located time-series (Fig.4a). A subsequent prompt 5 Figure 3: Validation scatter plot comparing in-situ temperature observations (0–500 m) from ten HAUSGARTEN moorings with co-located values from the Copernicus GLORYS12V1 reanaly- sis. The Oceanographer Agent retrieved dailythetao elds, performed 4D nearest-grid-point matchup (N= 135,678), and calculated statistics (Bias +0.35 ◦ C, RMSE 1.09 ◦ C,r= 0.31). Color scale indicates observation depth. directed the system to generate a wind rose diagram (Fig. 4b). The agents recovered from two runtime errors, including an out-of-bounds interpolationedge case (Supplementary Note 3). During interpolation of the ship track against the gridded ERA5 eld, the Oceanographer Agent detected that 83 of the 2,013 co-located points returned invalid values where the vessel trajectory grazed the bounding-box edge. Rather than reporting a partial result, the agent autonomously expanded the retrieval domain by approximately 0.3 ◦ and re-executed the extraction, achieving 100% spatial coverage. The wind rose generation required the agent to autonomously diagnose a dimensional mismatch between the number of speed-bin edges and the corresponding Beaufort scale labels (an o -by-one error inherent in binned classi cation) and correct the bin de nitions before producing the nal diagram. 2.2.4 Scenario 4: Biodiversity–water mass coupling Interdisciplinary research often requires merging sparsebiological observations with continu- ous physical elds. We tasked the system to compare jelly shbiodiversity between the Nor- wegian and Irminger Seas using a cruise dataset from the NorthAtlantic (May 2013; DOI: 10.1594/PANGAEA.829702) and to link species abundance to speci c water masses. First, the Ecologist Agent strati ed sampling stations into “Eastern” (Norwegian Sea) and “Western” (Irminger Sea) transects based on a user-de ned 20 ◦ W longitude threshold. It then computed the Shannon-Wiener diversity index (H ′ ) for each station. A Welch’s t-test (t= 6 (a) In-situ air temperature (blue) measured aboard RV Po- larstern compared with co-located ERA5 2-meter temperature (red) interpolated to the ship’s hourly positions. The system calculated the bounding box to minimize cloud data retrieval (∼260 MB). (b) Wind rose diagram for the PS122/1 cruise leg (September– December 2019,n= 120,616 observations). Wind directions are binned into 16 sectors; color scale indicates Beaufort wind force cate- gories. The diagram was produced after one self-correction cycle. Figure 4:ERA5/MOSAiC Lagrangian validation and wind regime characterization. −3.39,p= 0.0014) indicated signi cantly higher biodiversity in the Western transect ( ̄ H ′ west = 0.81) compared to the Eastern transect ( ̄ H ′ east = 0.35). During the diversity computation, the Ecologist Agent encountered a division-by-zero error caused by zero-abundance taxa. It autonomously applied a masking operation to exclude these entries before the logarithmic transformation. The Visualization Agent produced a comparative violin plot, utilizing the re exive quality control loop to detect and x label overlaps across ve iterations. Subsequently, to determine the environmental niche of the dominant hydrozoanAglan- tha digitale, the Oceanographer Agent queried the Copernicus Marine Service to retrieve co- located potential temperature (θ) and salinity (S) pro les. It performed a 4D nearest-neighbor matchup between the global reanalysis model and the local biological samples. The resulting Temperature-Salinity (T-S) diagram (Fig. 5) shows that the highest abundances ofA. digitale cluster within the cold, fresh end-member (θ <2 ◦ C,S≈34.8–35.0 psu), corresponding to Sub- polar Mode Water carried by the Irminger Current. This is consistent with the species’ known anity for subarctic water masses in the North Atlantic [ Pantiukhin et al.,2024]. The warmer, more saline Norwegian Atlantic Current stations (θ >5 ◦ C,S >35.1 psu) hosted markedly lower abundances (Supplementary Note 4). 2.3 Search benchmark To evaluate the agentic retrieval architecture, we conducted a benchmark using a curated dataset of 100 natural language queries spanning six geoscienti c domains: physical oceanogra- phy, marine geology, biogeochemistry, paleoclimatology,marine ecology, and atmospheric sci- ence (Supplementary Note 5). Queries were manually assembled to represent realistic research needs and strati ed into ve complexity categories: single-parameter lookups, multi-constraint ltering, spatiotemporal subsetting, cross-domain integration, and open-ended discovery. Each query’s retrieval results were evaluated using a structured rubric scored on a 1–10 scale across ve metrics: Precision, Recall, Parameter Coverage, Metadata Quality, and Result Relevance. Evaluation was performed programmatically using Gemini 3 Pro Preview [ Google DeepMind, 2025] as an automated judge, following a structured rubric with detailed scoring criteria for 7 Figure 5:Bio-physical coupling generated by the Oceanographer Agent.The system augmented the PANGAEA biological dataset by retrieving and co-locating 4D hydrographic data (θ/S) from the Copernicus Marine Service GLORYS12V1 reanalysis. The T-S diagram highlights the co-occurrence ofAglantha digitalewith distinct water masses (log 10 -abundance). each metric. All search architectures were driven by GPT-5.2[OpenAI,2026]. We compared three con gurations: (1) Baseline (Elasticsearch BM25, raw keyword matching), (2) Simple LLM (single-pass translation of the natural language queryinto a structured Boolean request), and (3) Agentic Search (our ReAct-based iterative agent, which decomposes user intent into multiple Boolean query permutations and re nes them acrosssuccessive retrieval cycles). The Agentic Search strategy increased the overall mean scoreto 8.14/10 compared to the Simple LLM’s 5.46/10 and the Baseline’s 2.87/10 (Fig. 6). The most signi cant gains occurred in Precision (M1: 8.53) and Parameter Coverage (M3: 8.99). Beyond aggregate statistics, qual- itative analysis of individual queries reveals distinct failure modes that the agentic architecture resolves (Supplementary Note 6). Temporal reasoning.When queried for “salinity pro les from the Weddell Sea during winter 2013,” both the Baseline and Simple LLM failed to resolve the hemispheric seasonal constraint: the former returned IceBird altimetry campaigns from the wrong year, while the latter selected summer cruises. Only the Agentic Search correctly mapped “Austral Winter” to June–August and ranked the ANT-XXIX/6 (AWECS) cruise as the top result (Agentic: 9.8/10 vs. Baseline: 1.8/10; full analysis in Supplementary Note 6). Source disambiguation.A query for “satellite altimetry data” exposed a systematicre- trieval error: both the Baseline and Simple LLM returned in-situ mooring datasets whose meta- data mentioned “used to validate satellite altimetry,” confusing ground-truth with the requested product. The Agentic Search correctly prioritized gridded sea-surface height anomaly products from CryoSat-2, Envisat, and TOPEX/Poseidon (Agentic: 9.6/10 vs. Baseline: 1.0/10). These examples indicate that iterative query re nement resolves implicit scienti c con- straints (hemispheric seasonality, measurement geometry, product-level intent) that keyword and single-pass approaches systematically miss. 8 M1: PrecisionM2: Spatio-TemporalM3: ParameterM4: DiversityM5: Ranking Evaluation Metric 0 2 4 6 8 10 12 Mean Score (1-10) 2.8 3.0 3.1 2.7 2.7 4.7 5.5 5.1 3.8 4.7 8.5 8.0 9.0 6.9 8.3 Search Method Performance Comparison (100 Queries) Elastic SearchSimple LLMAgentic Search Figure 6: Grouped bar chart comparing Baseline (Elastic Search), Simple LLM, and Agentic Search architectures across ve semantic metrics using 100geoscienti c queries. Error bars indicate standard deviation. The Agentic Search method outperforms traditional approaches in all categories. 3 Discussion Four architectural primitives together enable the composition of multi-step geoscienti c work- ows across heterogeneous data formats: data-type-aware routing, sandboxed execution, self- correction from execution tracebacks, and re exive visualquality control. Iterative agentic search treats user queries not as strings to be matched but asunder-speci ed scienti c intents to be progressively re ned. These primitives generalize across the validation scenarios without scenario-speci c engineering, and by resolving implicit constraints that keyword and single-pass LLM approaches miss, the system reduces the manual e ort required to transition from data retrieval to analysis [ Singh et al.,2025,Zhang et al.,2023a]. Agentic search works because it replicates what a domain expert does implicitly during manual retrieval. When a researcher searches for “sediment cores from the Fram Strait,” they do not submit these words verbatim; they translate the geographic name into a bounding-box speci cation, identify relevant campaign identi ers, andcross-reference parameter inventories. In practice, however, no user formulates this as a structured Boolean query against an Elas- ticsearch index; they type natural language and expect relevant results. The LLM bridges this gap by translating free-text intent into structured Boolean queries, and the agentic layer goes further by generating multiple query permutations anditeratively re ning them against the returned metadata. The four validation scenarios stress-test di erent capabilities of this approach: cross-domain data fusion (Scenario 1), format harmonization across heterogeneous tabular sources (Scenario 2), resource-ecient cloud-native retrieval along non-Eulerian trajec- tories (Scenario 3), and biological–physical coupling requiring external reanalysis integration (Scenario 4). We note that the system functions as an assisted analysis tool whose output quality depends on the speci city of the user’s intent; vague or under-constrained queries may yield retrieval of marginally relevant datasets. The layered self-correction strategy addresses error propagation at two levels. At the pro- grammatic level, tracebacks from the sandboxed Python kernel are fed back to the generating agent, which parses the error, refactors the code, and re-executes. In Scenario 1, for exam- ple, the Oceanographer Agent diagnosed and corrected a Copernicus API rejection caused by an invalid depth coordinate. At the perceptual level, the visual quality-control loop subjects every generated gure to a structured critique. In Scenario4, this loop drove the violin plot 9 through ve successive re nements (scores 3→4→4→5→9/10), each targeting a distinct ren- dering defect. Across all scenarios, the system resolved all ve runtime errors without human intervention, and the visual QC loop executed 11 critique cycles across the generated gures (Supplementary Table S6). This combination of programmatic and perceptual feedback draws on methods from verbal reinforcement learning [Du et al.,2024,Shinn et al.,2023] and multi- agent collaborative veri cation [ Darwish et al.,2025,Shi et al.,2025]. A third mechanism, cross-model escalation to a secondary foundation model (Section 5.2.5), provides an additional fallback, though it was not triggered in the four validationscenarios. However, the centralization of task decomposition in the Supervisor Agent constitutes a single point of failure. If the Supervisor incorrectly decomposes a query (for example, routing a gridded reanalysis comparison to the DataFrame Agent instead of the Oceanographer Agent), downstream agents have no mechanism to detect or correct this misassignment. Similarly, the system cannot verify domain-speci c correctness beyond code executability. For example, it cannot detect whether a co-located dataset contains a systematic spatial o set from an incorrect coordinate reference system, or whether an unnecessary unit conversion was applied. Addressing these failure modes will require domain-speci cvalidation hooks, such as automated CRS consistency checks and unit-conversion veri cation, that operate independently of the Supervisor’s planning logic. The Supervisor-Worker topology also facilitates reproducibility and extensibility. Each an- alytical step is captured as executable Python code within adeterministic sandbox, providing a complete provenance trail that complements the stochastic reasoning layer [ Herschel et al., 2017]. Extending the system to new domains requires three steps:de ning a domain-speci c agent prompt, registering routing rules that map feature ags to the new agent, and adding the required feature ags to the metadata introspection module (Section5.1.4). Natural can- didates include sedimentology (grain-size distributions, stratigraphic correlation), paleoclima- tology (proxy calibration, age–depth modelling), and geochemistry (elemental ratios, isotope systematics). This modular design follows principles established in prior multi-agent frame- works [Shen et al.,2023,Hong et al.,2024]. The current system operates reactively, responding to individual user queries. However, the architecture is compatible with an asynchronous mode in which agents traverse the archive over extended periods, systematically loading datasets, cross-referencing observations against reanal- ysis elds, and agging statistical anomalies, metadata inconsistencies, or previously unexam- ined inter-dataset correlations. For example, an agent could compare all CTD pro les in a given region against the corresponding reanalysis grid cells, identifying systematic observation–model discrepancies that may indicate instrument calibration drift or ocean features not captured by the model. A complementary application is inverse metadatageneration: agents could ingest poorly documented legacy les, infer parameter names and spatiotemporal coverage from data distributions and le headers, and produce machine-actionable metadata records, increasing the fraction of discoverable datasets [ Jacobsen et al.,2020]. Such long-running exploration would generate candidate ndings for subsequent human review, turning the archive from a passive store into one that actively surfaces candidate patterns for expert review. We note that this direction is speculative and would require advances in hallucination detection [ Farquhar et al., 2024], automated signi cance ltering to suppress spurious correlations, and human-in-the-loop validation protocols before operational deployment. Modularity is also a practical design requirement for systems of this kind. Because each worker agent communicates with the Supervisor only througha standardized message interface, individual agents can be detached from the pipeline and reused independently. The Oceanogra- pher Agent, for example, already integrates ERA5 and Copernicus data without relying on the PANGAEA retrieval layer, and could serve as a standalone reanalysis tool in other work ows. Designing agents to be self-contained from the outset lowers the cost of adapting the system to 10 new repositories or embedding individual components into third-party platforms. 4 Limitations and Future Directions Our results must be interpreted within the constraints of the underlying data infrastructure. The system’s retrieval performance is tied to the quality ofexisting metadata; in the PANGAEA archive, only approximately 15% of datasets contain structured abstracts, and sparse parameter descriptions limit the ecacy of semantic search for the remaining majority. This creates a “long-tail” problem where poorly documented datasets remain e ectively invisible to the retrieval layer. Additionally, while the architectural components are not speci c to PANGAEA, our evaluation was restricted to this repository. Transferability to other major repositories, such as NASA EOSDIS or GBIF, remains untested and will require handling distinct API conventions, authentication mechanisms, and metadata schemas. Because the system relies on stochastic inference, outputsare inherently non-deterministic, which con icts with the reproducibility requirements of scienti c computing. Repeated exe- cution of the same query may yield functionally equivalent but syntactically di erent code, di erent dataset retrieval orderings, and marginally di erent numerical results due to alterna- tive interpolation strategies or coordinate-matching heuristics. While the deterministic sandbox ensures that any given code block produces identical output, the upstream reasoning layer in- troduces variability in which code is generated. The system also remains vulnerable to hallucinations at twolevels [ Huang et al.,2023]. Visual hallucinations (rendering artefacts, incorrect axis labels, missing legends) are mitigated by the re exive quality-control loop described above. Interpretive hallucinations are more dif- cult to detect: the system might correctly compute a correlation coecient but attribute it to a causal mechanism unsupported by the data, or generate a narrative that over-interprets a statistically marginal result. The system also shows reduced performance on multi-step rea- soning that requires domain-speci c knowledge; for example, inferring that a depth coordinate mismatch implies an incompatible vertical reference framerather than a unit conversion error [ Wang et al.,2024]. As domain-speci c foundation models mature, integratingthem into the reasoning layer may reduce these failure modes. Each multi-step analytical session involves dozens of LLM inference calls spanning task planning, code generation, error diagnosis, and visual critique, resulting in non-trivial compu- tational costs when using commercial models such as GPT-5.2. Emerging open-weight models and model distillation techniques may reduce these costs, particularly for routine operations such as code refactoring and metadata parsing, enabling deployment on institutional computing clusters. From a security perspective, although the sandbox provides lesystem and process isolation, the system does not currently perform static analysis or formal veri cation of agent-generated code. Future hardening should include AST-level code inspection, per-session resource quotas (CPU time, memory, disk), and an explicit allowlist of importable Python packages [ Vangala et al. ,2025]. In multi-tenant deployments, session isolation must also prevent cross-user data leakage through shared kernel state. Finally, our evaluation methodology introduces a circularity: an LLM (Gemini 3 Pro Pre- view) serves as the automated judge for a system driven by another LLM (GPT-5.2). The judge evaluates factual alignment between returned datasets and query constraints but does not assess the reasoning process that led to dataset selection. While this approach enables reproducible scoring at scale, future validation must incorporate human expert annotation on a representative subset of queries to calibrate automated scores, quantify inter-rater agreement, and identify systematic blind spots in the automated evaluation. The full scoring rubric and per-query results are provided in Supplementary Note 5. 11 5 Methods 5.1 Search and retrieval To bridge the semantic gap between natural language inquiryand the PANGAEA archive schema, we designed and evaluated a multi-tiered retrievalarchitecture comprising three con- gurations of increasing autonomy (Supplementary Note 5). Each tier was benchmarked against a curated dataset of 100 natural language queries spanning six geoscienti c domains, scored by an automated judge on ve semantic metrics (Supplementary Note 5). 5.1.1 Baseline: Elasticsearch keyword retrieval The lowest tier relies on the PANGAEA Elasticsearch index usingstandard BM25 scoring with tf-idf weighting. Queries are matched verbatim against metadata elds (title, abstract, parameter names, campaign identi ers). This con gurationprovides sub-second latency but lacks semantic understanding: synonymous terms (e.g., “CTD” vs. “conductivity–temperature– depth”), geographic colloquialisms (“Fram Strait” vs. a bounding-box speci cation), and tem- poral abstractions (“Last Glacial Maximum”) are not resolved. The Baseline achieved a mean score of 2.87/10 in our benchmark, re ecting high recall fortrivially speci ed queries but poor precision and parameter coverage for complex scienti c requests. 5.1.2 Intermediate: simple LLM query translation The intermediate tier interposes a single LLM inference call between the user and the search engine. GPT-5.2, operating in a zero-shot setting, translates the natural language query into a structured Lucene Boolean string, resolving geographic names to coordinate ranges and expand- ing parameter synonyms. This single-pass approach improved the benchmark mean to 5.46/10 but remains brittle: the model has no mechanism to inspect the returned results, detect low- quality matches, or iteratively re ne the query. Ambiguous queries that require multiple rounds of disambiguation therefore saturate at moderate precision. 5.1.3 Advanced: agentic search The Search Agent employs a ReAct (Reasoning + Acting) loop [ Yao et al.,2023], operating through four stages: 1.Inferential Query Translation:The agent decomposes user intent into metadata con- straints, deriving implicit boundary conditions where needed; for example, translating “Fram Strait” or “Last Glacial Maximum” into bounding-box coordinates and numeric temporal intervals. 2.Execution and Introspection:The agent executes a structured query and evaluates the returned metadata. 3.Parallelized Query Expansion:Rather than relying on a single search string, the agent programmatically generates a matrix of Boolean querypermutations targeting dis- tinct metadata elds, permuting instrument types, parameter synonyms, and platform identi ers. These queries are executed concurrently against the archive. 4.Semantic Consolidation:The system aggregates results, performs deduplication based on unique identi ers, and ranks each result by scoring its metadata against prede ned relevance criteria derived from the original user prompt. 12 5.1.4 Metadata retrieval and feature introspection Before any analytical code is generated, the system pre-fetches the full metadata record for each selected dataset, including parameter inventories, matrix dimensions, coordinate refer- ence systems, and temporal extents. From these records, it derives binary feature ags that characterize each dataset along ve dimensions: whether it contains geospatial coordinates, whether it is stored as a gridded array, whether it includes adepth axis, whether its le size exceeds a con gurable threshold, and whether it carries a temporal index. The system injects these ags into the Supervisor’s planning context, where they directly determine the routing decision (Section 5.2.2). For instance, a dataset agged as gridded with a depth axisis routed to the Oceanographer Agent, whereas a at table without spatial coordinates is assigned to the DataFrame Agent. This pre-execution introspection eliminates the ambiguity that would otherwise arise if the Supervisor inferred data type from the user’s natural language description alone. 5.2 Multi-agent orchestration The analytical core is a stateful MAS governed by a directed cyclic graph, implemented using the LangGraph framework [ LangChain,2024b] built on top of LangChain [LangChain,2024a], with a Supervisor Agent as the orchestrator. The architecture isolates the reasoning engine from the computational kernel, enabling task routing basedon data-type introspection and ensuring reproducibility through deterministic code execution. 5.2.1 Supervisor agent and task planning The Supervisor maintains the global conversational state and manages inter-agent handovers. Before delegating tasks, it inspects the internal structure of retrieved resources (e.g., NetCDF dimension hierarchies, Zarr chunk layouts) to select the appropriate specialist. It decomposes user queries into sub-tasks via a planning module, dispatches them sequentially, and updates the plan as results become available. Safeguards include aniteration cap and strict data- type routing constraints to prevent, for instance, routingmulti-dimensional arrays to agents specialized for tabular data. 5.2.2 Data-type routing and agent specialization The introspection performed during the metadata enrichment phase (Section 5.1.4) determines the execution branch, routing tasks to speci c worker nodesbased on their architectural spe- cialization: •The Oceanographer Agent:Specialized forN-dimensional array handling (NetCDF/Zarr). It possesses tool-augmented access to external climate reanalysis APIs (Copernicus Marine Service, ERA5 via Arraylake [ Earthmover,2025]), enabling contextualization of local observations within global environmental elds. The agent automatically extracts spatiotemporal bounds from loadeddatasets, retrieves the match- ing reanalysis subset via lazy-loaded cloud-native Zarr stores, and performs geometric co-location (nearest-grid-point for Eulerian comparisons, 4D interpolation for Lagrangian trajectories). Full details are provided in SupplementaryNote 7. •The Ecologist Agent:Optimized for biodiversity metrics and tabular data manipula- tion, equipped with specialized libraries for biological indices calculation. 13 •The Visualization Agent:A general-purpose plotting specialist augmented with a Retrieval-Augmented Generation (RAG) index of domain-speci c exemplar plots, en- abling it to match gure types to data dimensionality and enforce scienti c conventions (e.g., reversed depth axes, standard color palettes). •The DataFrame Agent:Strictly con ned to data manipulation, statistical analysis, and ltering of tabular data. This agent is architecturallyprohibited from generating visualizations, enforcing a clean separation between computation and presentation. •The Writer Agent:A synthesis specialist that collects outputs from upstreamagents and produces a coherent scienti c narrative, interpretinggenerated gures via a vision- language model interface. 5.2.3 Adaptive context management Building on the multi-tier memory concept described in Pantiukhin et al.[2025], we implement a progressive summarization strategy to manage context window limits. The system partitions conversational history into two segments: recent messagesretained at full resolution, and older exchanges compressed into a structured summary by a fast inference model. This summary explicitly preserves critical operational state (variable names, loaded dataset identi ers, le paths, and intermediate results) while discarding conversational ller. The result is that agents retain sucient context to resume multi-step work ows evenwhen the raw token history exceeds the model’s context window. 5.2.4 Sandboxed execution Agents synthesize Python code executed in a persistent Jupyter kernel within an isolated lesys- tem sandbox. This design ensures three critical properties. First,State Persistence: unlike stateless function calls, the kernel maintains active memory of variables, loaded datasets, and imported libraries, allowing iterative multi-step analysis without data reloading. Second,Ses- sion Isolation: each user session initializes a unique workspace identi ed by a universally unique identi er (UUID), with all data retrieval and artifact generation con ned to this directory to prevent cross-session contamination. Third,Self-Correction: the kernel captures standard error outputs (tracebacks) and feeds them back to the generating agent, which parses the traceback, refactors the code, and re-executes. 5.2.5 Cross-model error escalation (Wise Agent) When a worker agent fails to resolve an error after a con gurable number of retry cycles, the system escalates the task to a secondary agent (theWise Agent) backed by a di erent foundation model. In the current deployment, the primary pipeline runs on GPT-5.2 [ OpenAI, 2026], while the Wise Agent runs on Claude Opus 4.6 [Anthropic,2026]. The escalation transfers the full execution context (code, traceback, and prior repair attempts) to the Wise Agent, which generates an independent x. This cross-model fallback exploits the observation that di erent LLMs exhibit complementary failure modes: an error patternthat causes one model to loop often lies outside the failure distribution of another. In practice, this mechanism resolved deadlocked error cycles that single-model retries could not break. 5.2.6 Visual quality control Visual outputs are subjected to a re exive quality-control loop that intercepts every generated gure before it is surfaced to the user. Once a plot is saved tothe sandbox, the image is 14 passed to a vision-enabled language model acting as a structured critic. The critic evaluates the gure against a ten-point rubric spanning ve dimensions: axis labelling and units, legend completeness, adherence to domain conventions (e.g., inverted depth axes for oceanographic pro les, standard color palettes for temperature elds), absence of overlapping text, and overall readability at publication scale. If the composite score falls below a strict threshold, the critic returns structured, actionable feedback specifying each de ciency. This feedback is appended to the generating agent’s context, triggering a targeted code refactoring cycle. The loop repeats until the score exceeds the threshold or a maximum of ve iterations is reached. Data Availability All datasets analyzed in this study are publicly available through the PANGAEA data repository. Speci c datasets are referenced by DOI in the corresponding scenario descrip- tions: Weddell Sea microplastics ( 10.1594/PANGAEA.941076), MOSAiC PS122/1 navigation ( 10.1594/PANGAEA.935221), North Atlantic jelly sh survey (10.1594/PANGAEA.829702), and HAUSGARTEN moorings (PANGAEA DOIs 845608–845617). Climate reanalysis data were accessed from the Copernicus Marine Service (GLORYS12V1) and the ERA5 reanalysis accessed via the Arraylake cloud-native data platform (Earthmover, https://earthmover.io). The benchmark query dataset (100 queries, Supplementary Note 5) is available at https: //github.com/CliDyn/pangaeaGPT . Code Availability The PANGAEA-GPT source code, agent con gurations, and benchmark evaluation scripts are available at https://github.com/CliDyn/pangaeaGPT. Author Contributions D.P. conceived and designed the system, developed the multi-agent architecture, conducted all experiments, and wrote the manuscript. I.K. provided domain expertise in physical oceanog- raphy, guided the design of the validation scenarios, and reviewed the scienti c interpretation of the results. B.S. contributed to the system design and assisted with the integration of the PANGAEA Elasticsearch API and metadata harvesting protocols. A.A.J. contributed to the ecological use-case design, provided expertise in biodiversity metrics, and assisted with the curation of the benchmark query dataset. N.K. supervised theproject, provided strategic guid- ance on the system architecture and its integration with thePANGAEA infrastructure, and contributed to manuscript revision. T.J. provided institutional support, strategic oversight, and contributed to manuscript revision. All authors reviewed and approved the nal manuscript. Competing Interests The authors declare no competing interests. Acknowledgements This work was supported by the Helmholtz Association and the Federal Ministry of Educa- tion and Research (BMBF) for supporting the DataHub Initiative of the Research Field Earth 15 and Environment; the European Union’s Destination Earth Initiative and relates to tasks en- trusted by the European Union to the European Centre for Medium-Range Weather Forecasts implementing part of this Initiative with funding by the European Union; and Projects S1: Di- agnosis and Metrics in Climate Models of the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean,” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project no. 274762653). Views and opinions expressed are those of the authors only and do not necessarily re ect those of theEuropean Union or the European Climate Infrastructure and Environment Executive Agency (CINEA). Neither the European Union nor the granting authority can be held responsible for them. References Anthropic. Introducing claude opus 4.6. https://w.anthropic.com/news/claude-opus-4 -6 , 2026. Accessed: 2026-02-24. Zhen Bi, Ningyu Zhang, Yida Xue, Yixin Ou, Daxiong Ji, Guozhou Zheng, and Huajun Chen. Oceangpt: A large language model for ocean science tasks. InProceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 3325–3346, 2024. doi: 10.18653/v1/2024.acl-long.184. Daniil A Boiko, Robert MacKnight, Gabe Kline, and Gabe Gomes. Emergent autonomous scienti c research capabilities of large language models.arXiv preprint arXiv:2304.05332, 2023a. doi: 10.48550/arXiv.2304.05332. Daniil A Boiko, Robert MacKnight, Mark Kline, and Gabe Gomes. Autonomous chemical research with large language models.Nature, 624:570–578, 2023b. doi: 10.1038/s41586-023 -06792-0. Andres M Bran, Sam Cox, Oliver Schilter, Carlo Baldassari, Andrew D White, and Philippe Schwaller. Augmenting large-language models with chemistry tools.Nature Machine Intel- ligence, 6:525–535, 2024. doi: 10.1038/s42256-024-00832-8. Y Chen, W Wang, S Lobry, and C Kurtz. An llm agent for automatic geospatial data analysis. arXiv preprint arXiv:2410.18792, 2024. doi: 10.48550/arXiv.2410.18792. Ahmed M. Darwish, Essam A. Rashed, and Ghada Khoriba. Mitigating LLM hallucinations using a multi-agent framework.Information, 16(7):517, 2025. doi: 10.3390/info16070517. Yilun Du, Shuang Li, Antonio Torralba, Joshua B Tenenbaum, and Igor Mordatch. Improving factuality and reasoning in language models through multiagent debate. InInternational Conference on Machine Learning (ICML), 2024. doi: 10.48550/arXiv.2305.14325. Earthmover. Arraylake: Cloud-native data management for n-dimensional arrays. https: //earthmover.io/arraylake , 2025. Accessed: 2026-02-08. Sebastian Farquhar, Jannik Kossen, Lorenz Kuhn, and Yarin Gal. Detecting hallucinations in large language models using semantic entropy.Nature, 630:625–630, 2024. doi: 10.1038/s4 1586-024-07421-0. Janine Felden, Lars M”oller, Uwe Schindler, et al. Pangaea-data publisher for earth & envi- ronmental science.Scienti c Data, 10(1):347, 2023. doi: 10.1038/s41597-023-02269-x. 16 Alireza Ghafarollahi and Markus J Buehler. Atomagents: Alloydesign and discov- ery through physics-aware multi-modal multi-agent arti cial intelligence.arXiv preprint arXiv:2407.10022, 2024a. doi: 10.48550/arXiv.2407.10022. Alireza Ghafarollahi and Markus J Buehler. Sciagents: Automating scienti c discovery through multi-agent intelligent graph reasoning.arXiv preprint arXiv:2409.05556, 2024b. doi: 10.4 8550/arXiv.2409.05556. Ali Essam Ghareeb, Benjamin Chang, Ludovico Mitchener, et al. Robin: A multi-agent system for automating scienti c discovery.arXiv preprint arXiv:2505.13400, 2025. doi: 10.48550/a rXiv.2505.13400. Yolanda Gil, C’edric H David, Ibrahim Demir, Bakinam T Essawy, Robinson W Fulweiler, Jonathan L Goodall, Leif Karlstrom, Huikyo Lee, Heath J Mills, Ji-Hyun Oh, et al. Toward the geoscience paper of the future: Best practices for documenting and sharing research from data to software to provenance.Earth and Space Science, 3(10):388–415, 2016. doi: 10.1002/2015EA000136. Google DeepMind. Gemini 3. https://blog.google/products-and-platforms/products/ gemini/gemini-3/ , 2025. Accessed: 2026-02-08. Juraj Gottweis, Wei-Hung Weng, Alexander Daryin, et al. Towards an AI co-scientist.arXiv preprint arXiv:2502.18864, 2025. doi: 10.48550/arXiv.2502.18864. Siyuan Guo, Cheng Deng, Ying Wen, Hechang Chen, Yi Chang, and Jun Wang. DS-Agent: Automated data science by empowering large language models with case-based reasoning. In International Conference on Machine Learning (ICML), 2024a. doi: 10.48550/arXiv.2402. 17453. Taicheng Guo, Xiuying Chen, Yaqi Wang, et al. Large language model based multi-agents: A survey of progress and challenges. InProceedings of the Thirty-Third International Joint Conference on Arti cial Intelligence (IJCAI-24), pages 8039–8047, 2024b. doi: 10.24963/ijc ai.2024/890. Melanie Herschel, Ralf Diestelk ̈amper, and Houssem Ben Lahmar. A survey on provenance: What for? what form? what from?The VLDB Journal, 26(6):881–906, 2017. doi: 10.1007/ s00778-017-0486-1. Sirui Hong, Xiawu Zheng, Jonathan Chen, et al. Metagpt: Meta programming for multi-agent collaborative framework. InThe Twelfth International Conference on Learning Representa- tions (ICLR), 2024. doi: 10.48550/arXiv.2308.00352. Sirui Hong, Yizhang Lin, Bang Liu, Bangbang Liu, Binhao Wu, Danyang Li, Jiaqi Chen, Jiayi Zhang, Jinlin Wang, Li Li, et al. Data interpreter: An LLM agentfor data science. In Findings of the Association for Computational Linguistics: ACL 2025, 2025. doi: 10.48550 /arXiv.2402.18679. Lei Huang, Weijiang Yu, Weitao Ma, Weihong Zhong, Zhangyin Feng, Haotian Wang, Qiang- long Chen, Weihua Peng, Xiaocheng Feng, Bing Qin, and Ting Liu. A survey on hallucina- tion in large language models: Principles, taxonomy, challenges, and open questions.arXiv preprint arXiv:2311.05232, 2023. doi: 10.48550/arXiv.2311.05232. Annika Jacobsen, Ricardo de Miranda Azevedo, Nick Juty, et al. Fair principles: Interpretations and implementation considerations.Data Intelligence, 2(1-2):10–29, 2020. doi: 10.1162/dint r00024. 17 Ben Kobler, John Berbert, P Caulk, and P C Hariharan. Architecture and design of storage and data management for the nasa earth observing system data andinformation system (eosdis). InProceedings of IEEE 14th Symposium on Mass Storage Systems, pages 65–76. IEEE, 1995. doi: 10.1109/MASS.1995.528217. Nikolay Koldunov and Thomas Jung. Local climate services for all, courtesy of large language models.Communications Earth & Environment, 5:13, 2024. doi: 10.1038/s43247-023-01199 -1. Ivan Kuznetsov, Antonia Anna Jost, Dmitrii Pantiukhin, Boris Shapkin, Thomas Jung, and Nikolay Koldunov. Transforming climate services with LLMs and multi-source data integra- tion.npj Climate Action, 4:97, 2025. doi: 10.1038/s44168-025-00300-y. LangChain. Langchain: Build context-aware reasoning applications. https://github.com/l angchain-ai/langchain , 2024a. Accessed: 2026-02-24. LangChain. Langgraph: Build resilient language agents as graphs. https://github.com/lan gchain-ai/langgraph , 2024b. Accessed: 2026-02-24. Jean-Michel Lellouche, Eric Greiner, Romain Bourdall ́e-Badie, Gilles Garric, Ang ́elique Melet, Marie Dr ́evillon, Cl ́ement Bricaud, Mathieu Hamon, OlivierLe Galloudec, Charly Regnier, et al. The Copernicus global 1/12 ◦ oceanic and sea ice GLORYS12 reanalysis.Frontiers in Earth Science, 9:698876, 2021. doi: 10.3389/feart.2021.698876. Guohao Li, Hasan Abed Al Kader Hammoud, Hani Itani, Dmitrii Khizbullin, and Bernard Ghanem. CAMEL: Communicative agents for “mind” explorationof large language model society. InAdvances in Neural Information Processing Systems (NeurIPS), 2023a. doi: 10.48550/arXiv.2303.17760. Xin Li, Min Feng, Youhua Ran, et al. Big data in earth system science and progress towards a digital twin of earth.Nature Reviews Earth & Environment, 4:319–332, 2023b. doi: 10.103 8/s43017-023-00409-w. Zhouhan Lin, Cheng Deng, Le Zhou, et al. Geogalactica: A scienti c large language model in geoscience.arXiv preprint arXiv:2401.00434, 2023. doi: 10.48550/arXiv.2401.00434. Luke Harry Marsden, P ̊al Gunnar Ellingsen, Lara Ferrighi, etal. Best practices for data management in marine science: lessons from the Nansen Legacyproject.Earth System Science Data, 17:5983–5996, 2025. doi: 10.5194/essd-17-5983-2025. H Ning, Z Li, T Akinboyewa, and M N Lessani. An autonomous gis agent framework for geospatial data retrieval.International Journal of Digital Earth, 18:2458688, 2025. doi: 10.1080/17538947.2024.2458688. Alexander Novikov, Ngˆan Vu, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, et al. Alphaevolve: A coding agent for scienti c and algorithmic discovery.arXiv preprint arXiv:2506.13131, 2025. doi: 10.48550/a rXiv.2506.13131. Google DeepMind. OpenAI. Introducing gpt-5.2. https://openai.com/index/introducing-gpt-5-2/, 2026. Accessed: 2026-02-08. Dmitrii Pantiukhin, Gerlien Verhaegen, and Charlotte Havermans. Pan-Arctic distribution modeling reveals climate-change-driven poleward shifts of major gelatinous zooplankton species.Limnology and Oceanography, 69(6):1316–1334, 2024. doi: 10.1002/lno.12568. 18 Dmitrii Pantiukhin, Boris Shapkin, Ivan Kuznetsov, AntoniaAnna Jost, and Nikolay Koldunov. Accelerating earth science discovery via multi-agent LLM systems.Frontiers in Arti cial Intelligence, 8:1674927, 2025. doi: 10.3389/frai.2025.1674927. David Rein, Betty Li Hou, Asa Cooper Stickland, Jackson Petty, Richard Yuanzhe Pang, Julien Dirani, Julian Michael, and Samuel R Bowman. GPQA: A graduate-level google-proof q&a benchmark. InFirst Conference on Language Modeling, 2024. doi: 10.48550/arXiv.2311.12 022. Nicol’as Robinson-Garc’ıa, Evaristo Jim’enez-Contreras, and Daniel Torres-Salinas. Analyz- ing data citation practices using the data citation index.Journal of the Association for Information Science and Technology, 67(12):2964–2975, 2016. doi: 10.1002/asi.23529. Timo Schick, Jane Dwivedi-Yu, Roberto Dess`ı, Roberta Raileanu, Maria Lomeli, Eric Hambro, Luke Zettlemoyer, Nicola Cancedda, and Thomas Scialom. Toolformer: Language models can teach themselves to use tools. InAdvances in Neural Information Processing Systems (NeurIPS), volume 36, pages 68539–68551, 2023. doi: 10.48550/arXiv.2302.04761. John L Schnase, Daniel Q Du y, Glenn S Tamkin, Denis Nadeau, JohnH Lindholm, Michael A McInerney, and WP Webster. Big data challenges in climate science: Improving the next- generation cyberinfrastructure.IEEE Geoscience and Remote Sensing Magazine, 4(3):10–22, 2016. doi: 10.1109/MGRS.2016.2574107. Yongliang Shen, Kaitao Song, Xu Tan, Dongsheng Li, Weiming Lu, and Yueting Zhuang. HuggingGPT: Solving AI tasks with ChatGPT and its friends in Hugging Face. InAdvances in Neural Information Processing Systems (NeurIPS), 2023. doi: 10.48550/arXiv.2303.17580. Jinxin Shi, Jiabao Zhao, Xingjiao Wu, Ruyi Xu, Yuan-Hao Jiang, and Liang He. Mitigat- ing reasoning hallucination through multi-agent collaborative ltering.Expert Systems with Applications, 263:125723, 2025. doi: 10.1016/j.eswa.2024.125723. Noah Shinn, Federico Cassano, Ashwin Gopinath, Karthik Narasimhan, and Shunyu Yao. Re- exion: Language agents with verbal reinforcement learning. InAdvances in Neural Infor- mation Processing Systems (NeurIPS), 2023. doi: 10.48550/arXiv.2303.11366. Aditi Singh, Abul Ehtesham, Saket Kumar, and Tala Khoramshahi. Agentic retrieval- augmented generation: A survey on agentic RAG.arXiv preprint arXiv:2501.09136, 2025. doi: 10.48550/arXiv.2501.09136. Bhanu Prakash Vangala, Ali Adibifar, Ashish Gehani, and Tanu Malik. AI-generated code is not reproducible (yet): An empirical study of dependency gaps in LLM-based coding agents. arXiv preprint arXiv:2512.22387, 2025. doi: 10.48550/arXiv.2512.22387. Xiaoxuan Wang, Ziniu Hu, Pan Lu, Yanqiao Zhu, Jieyu Zhang, Satyen Subramaniam, Arjun R Loomba, Shichang Zhang, Yizhou Sun, and Wei Wang. SciBench: Evaluating college-level scienti c problem-solving abilities of large language models. InInternational Conference on Machine Learning (ICML), 2024. doi: 10.48550/arXiv.2307.10635. Mark D Wilkinson, Michel Dumontier, IJsbrand Jan Aalbersberg, et al. The fair guiding principles for scienti c data management and stewardship.Scienti c Data, 3:160018, 2016. doi: 10.1038/sdata.2016.18. Yutaro Yamada, Varun Nagarajan, Jiaxin Song, et al. The ai scientist-v2: Workshop-level automated scienti c discovery via agentic tree search.arXiv preprint arXiv:2504.08066, 2025. doi: 10.48550/arXiv.2504.08066. 19 Shunyu Yao, Je rey Zhao, Dian Yu, Nan Du, Izhak Shafran, Karthik Narasimhan, and Yuan Cao. React: Synergizing reasoning and acting in language models.arXiv preprint arXiv:2210.03629, 2023. doi: 10.48550/arXiv.2210.03629. Wenqi Zhang, Yongliang Miao, Yueqi Li, Pei Deng, Dongmei Zhang, Tao Lv, Yaobo Li, Jianfeng Gao, Li Dong, and Furu Wei. Data-copilot: Bridging billionsof data and humans with autonomous work ow.arXiv preprint arXiv:2306.07209, 2023a. doi: 10.48550/arXiv.2306. 07209. Yifan Zhang, Cheng Wei, Shangyou Wu, Zhengting He, and Wenhao Yu. Geogpt: Un- derstanding and processing geospatial tasks through an autonomous gpt.arXiv preprint arXiv:2307.07930, 2023b. doi: 10.48550/arXiv.2307.07930. 20 Supplementary Material A Hierarchical Multi-Agent System for Autonomous Discovery inGeoscienti c Data Archives Supplementary Note 1: Scenario 1 extended protocol This section provides the complete methodological detail for Scenario 1 (Weddell Sea microplas- tic exploration), including verbatim user prompts, the agent task decomposition, autonomous error-recovery episodes, and the visualization re nementcycle. The condensed main-text ac- count (Scenario 1 in the main text) should be read alongside Supplementary Table S1 below. 1.1 Data and scienti c context The user loaded a published dataset of sea-surface microplastic concentrations collected by Manta trawl during two RV Polarstern expeditions to the Weddell Sea (PS111, austral summer 2018; PS117, austral summer 2019; Dataset DOI: 10.1594/PANGAEA.941076). The dataset comprisesN= 34 sampling stations spanning 78 ◦ S–60 ◦ S and 57 ◦ W–12 ◦ E, with microplastic abundance reported as particle counts per unit area (items km −2 ). 1.2 User prompts and task decomposition The work ow was driven by three user prompts. The rst directed the system to retrieve Copernicus Marine surface currents, calculate mean current speed over austral summers 2018– 2019, and create a regional map. The Supervisor decomposed this into two sub-tasks: (i) data retrieval and temporal averaging, routed to the Oceanographer Agent; and (i) cartographic rendering, routed to the Visualization Agent. The second prompt requested a minor aesthetic re nement (40% transparency). The third prompt directed the statistical analysis (Spearman correlation, Mann–WhitneyUtest, log-linear regression plot). 1.3 Self-correction episodes Two autonomous error-recovery events occurred during execution: Episode 1: Copernicus depth coordinate error (steps 4–5).The rst API request speci edminimum depth=0.0, which the Copernicus Marine Service rejected because the depth axis begins at 0.494 m. The agent parsed the error message, extracted the minimum valid depth coordinate, corrected the parameter, and resubmitted bothrequests successfully. Episode 2: Visualization layout failure (step 8).The re exive quality control loop rated the initial regional map at 3/10. The Visualization Agent autonomously refactored the gure layout, achieving a revised score of 9/10. 1 1.4 Statistical analysis protocol The DataFrame Agent co-located each of the 34 sampling stations with the nearest grid point of the mean current speed eld, producing a matched dataset ofN= 21 stations after ltering. Statistics: Spearmanρ=−0.466 (p= 0.033); Mann–WhitneyU= 88.0 (p= 0.022); OLS regression on log 10 (abundance): slope =−4.03, intercept = 3.58. Table S1: Annotated execution log for Scenario 1 (Weddell Sea microplastic exploration). Bold text indicates autonomous agent decisions. The work ow comprises three user prompts (Steps 1, 9, 11) and eleven autonomous agent actions, including one API error-recovery episode (Steps 4–5) and one re exive visualization re nement cycle (Step 8). # ActorActionContent / Reasoning 1UserPrompt“Retrieve Copernicus Marine surface currents (uo,vo) for [78 ◦ S–60 ◦ S, 57 ◦ W–12 ◦ E]. Calculate mean current speed aver- aged over austral summers 2018–2019. Create a regional map with current speed heatmap, streamlines, and sampling stations sized by microplastic abundance.” 2SupervisorPlanningDecomposes into two sub-tasks: (i) data retrieval and process- ing→Oceanographer Agent; (i) cartographic visualization→ Visualization Agent. 3OceanographerTemporal infer- ence Autonomouslyparses theDate/Timecolumn of the obser- vational dataset and determines the sampling campaigns span January–February 2018 and January–February 2019. Con- structs two separate API requests with matching temporal bounds. 4Oceanographer API call (fail)First Copernicus Marine request speci esminimumdepth=0.0. Service rejects the request: depth axis begins at 0.494 m. 5OceanographerSelf-correctionParses the error message, extracts the minimum valid depth co- ordinate (0.494 m), and resubmits. Both requests (2018, 2019) succeed. Downloads dailyu,v elds for the speci ed domain. 6Oceanographer Code executionComputes temporal mean of daily velocity elds across both summers usingxarray. Derives current speed as √ u 2 +v 2 . Saves result as NetCDF (mean surfacecurrentsJF2018- 2019 Weddell.nc). 7Visualization Code executionGenerates map using South Polar Stereographic projection (cartopy). Overlays streamlines on current speed heatmap (cmocean.deepcolormap). Plots stations as circles scaled by abundance, color-coded by concentration (magmacolormap). 8VisualizationRe ection loopRe exive quality control scores initial gure at3/10: “leg- end overlaps with map area, signi cantly impacts readability.” Agent refactors layout: repositions colorbars into dedicated axes, relocates size legend outside the plot area, adjusts sub- plot margins. Revised version scores9/10: “no overlapping elements, e ective data communication.” 9UserRe nement“Make the same gure, but draw circles with microplastics that are 40% transparent.” 10 Visualization Code executionRe-executes the plotting script with a single parameter change (alpha=0.6). 11 UserAnalysis request“Calculate the Spearman Rank Correlation betweencurrent speed and microplastic abundance. Perform a Mann–Whitney Utest splitting stations into Low vs. High Kinetic Energy groups at the median speed. Generate a log-linear regression plot with 95% CI and annotated statistics.” 2 # ActorActionContent / Reasoning 12 DataFrameCo-locationLoads the mean current speed eld. For each of the 34 sam- pling stations, extracts the nearest-grid-point current speed via xarrayvectorized selection. Filters toN= 21 stations with positive abundance and valid current speed. 13 DataFrameStatisticsComputes Spearmanρ=−0.466 (p= 0.033). Splits sta- tions at median speed (0.034 m s −1 ):n low = 11,n high = 10. Mann–WhitneyU= 88.0 (p= 0.022). Fits OLS regression on log 10 (abundance): slope =−4.03, intercept = 3.58. 14 Visualization Code executionGenerates scatter plot with logarithmicy-axis, regression line, 95% CI band, and statistical annotations. Re exive quality con- trol scores 8/10 on rst attempt. Supplementary Note 2: Scenario 2 extended protocol This section documents the complete interaction sequence for Scenario 2 (Copernicus/HAUSGARTEN mooring validation), including data characteristics, co-location methodology, and all au- tonomous error-recovery episodes. 2.1 Data and scienti c context The user loaded ten published mooring temperature datasetsfrom the LTER HAUSGARTEN observatory in the eastern Fram Strait (PANGAEA DOIs 10.1594/PANGAEA.845608 through 10.1594/PANGAEA.845617). Together the datasets contain 340,022 records spanning 2006– 2012, with instrument depths from 81–2,737 m. 2.2 Self-correction episodes One autonomous error-recovery event occurred: Episode 1: Datetime type mismatch (Oceanographer agent).AttributeError: numpy.datetime64has no.date()method. The agent refactored topd.Timestamp.normalize() for consistent datetime keys. 2.3 Validation statistics Overall (N= 135,678): Bias +0.35 ◦ C, RMSE 1.09 ◦ C, Pearsonr= 0.31. Depth binNBias ( ◦ C) RMSE ( ◦ C)r 50–100 m 7,738+0.200.97−0.05 100–200 m 47,119+0.421.06+0.11 200–500 m 80,821+0.321.12+0.42 Table S2: Depth-binned validation statistics for GLORYS12V1potential temperature against HAUSGARTEN mooring observations (0–500 m, 2006–2012). 3 Table S3: Annotated execution log for Scenario 2 (HAUSGARTEN mooring validation). Bold text indicates autonomous agent deci- sions. The work ow comprises three user prompts (Steps 1, 9, 11) and sixteen autonomous agent actions, including one error-recovery episode (Steps 15–16, datetime type mismatch) and two re exive visualization quality control cycles (Steps 6–8, 19). # ActorActionContent / Reasoning 1UserPrompt“Generate a plot visualizing the distribution of instrument depths across all loaded mooring datasets.” 2SupervisorPlanningDecomposes into two sub-tasks: (i) multi-dataset ingestion and depth extraction→Visualization Agent (with inline data wran- gling); (i) exploratory depth-distribution gure→Visualiza- tion Agent. 3VisualizationData pro lingLoads all 10 CSV les and inspects column names. Identi es Depth wateras the depth column across all datasets. Notes that three datasets (FEVI21, FEVI22, FEVI23) contain an ad- ditionalSalcolumn (14 columns vs. 13). Records total row counts: 340,022 records across 10 moorings. 4VisualizationData engineer- ing Computes per-instrument unique depths to avoid time-repeated values: groups byGear IDwithin each event and takes the median ofDepth water. Produces a consolidated table ofN= 35 unique instrument depths across 10 moorings, spanning 81– 2,737 m. 5Visualization Code executionGenerates a two-panel gure: (top) overall histogram with KDE of the 35 unique instrument depths; (bottom) violin + box + strip plot faceted by mooring event, ordered by median depth. Annotates per-mooring instrument counts (n= 3–5). 6VisualizationRe ection loopRe exive quality control scores the initial gure at7/10: “Axis labels slightly small; histogram x-axis extends into negative val- ues; font sizes could be increased for readability.” 7VisualizationSelf-correctionRefactors the gure: increases font sizes to 14 pt for axis la- bels and 16 pt for titles; clips KDE at zero (clip=(0, None)); setsxlim(left=0); xes seabornFutureWarningby explicitly assigninghueparameter withlegend=False. 8VisualizationRe ection loopRevised gure scores9/10: “All labels clear, no overlaps, sci- enti c conventions followed, depth represented correctly on the Y-axis.” Final plot accepted. 9UserPrompt“Create a consolidated dataset by merging all loaded mooring records from the upper ocean, speci cally between 0 and 500 meters depth. Preserve all original data columns and save to results.” 10 DataFrameCode executionIterates over all 10 datasets: coercesDepth waterto numeric, lters to 0≤Depth water≤500 m, appends amooring idcol- umn, and concatenates. Produces merged table: 136,183 rows× 15 columns. Savesmerged mooringsupperocean0500m.csv. 11 UserPrompt“Validate Copernicus global reanalysis data against in-situ mooring observations ( rst 500 m). Determine spatio-temporal bounds, retrieve dailythetao, co-locate with moorings, calcu- late Bias/RMSE/Correlation, and generate a scatter plot. Do it in one step with the Oceanographer Agent.” 4 # ActorActionContent / Reasoning 12 Oceanographer Code executionLoads all 10 mooring datasets, merges (340,022 rows), parses Date/Timeto UTC datetime, coerces numeric columns. Filters to 0–500 m with valid temperature, time, latitude, and longi- tude:N= 135,678 records. Extracts bounds: time 2006-09-05 to 2012-07-23; lat 78.61 ◦ N–79.74 ◦ N; lon 4.33 ◦ E–5.16 ◦ E; depth 81–500 m. 13 Oceanographer API callSubmits Copernicus Marine Service request forthetaofrom GLORYS12V1 (cmems modglophymy0.083degP1D-m) with the extracted spatio-temporal bounds. Receives NetCDF: 2,149 daily time steps×31 depth levels×13 latitude×10 longitude grid points. 14 Oceanographer Code executionImplements 4D nearest-grid-point co-location: builds a date→index lookup dictionary from the model time axis; com- putes nearest-neighbor indices for latitude, longitude, and depth using vectorized operations. Extracts modelthetaovalues at all 135,678 observation points via NumPy advanced indexing. 15 Oceanographer Code execution (fail) AttributeError:numpy.datetime64object has no at- tribute.date()during time-index mapping. The agent’s dictionary comprehension attempted to call.date() onnumpy.datetime64elements returned bypd.to - datetime().values. 16 OceanographerSelf-correctionIdenti es the type mismatch between NumPy date- time64 and Python native datetime. Refactors to use pd.Timestamp.normalize()for both the model time axis and the observation dates, producing consistentTimestamp keys. Re-executes co-location: 100% valid fraction (all 135,678 matchups successful). 17 Oceanographer StatisticsComputes overall validation statistics: Bias = +0.349 ◦ C, RMSE = 1.093 ◦ C, Pearsonr= 0.310. Computes depth- bin statistics: 50–100 m (n= 7,738, Bias +0.20, RMSE 0.97, r=−0.05); 100–200 m (n= 47,119, Bias +0.42, RMSE 1.06, r= +0.11); 200–500 m (n= 80,821, Bias +0.32, RMSE 1.12, r= +0.42). 18 Oceanographer Code executionGenerates validation scatter plot: mooring temperature vs. GLORYS12V1thetao, color-coded by observation depth (viridiscolormap), with 1:1 reference line and annotation box displayingN, Bias, RMSE,r. Subsamples to 50,000 points for rendering. 19 OceanographerRe ection loopRe exive quality control scores the scatter plot at8/10: “Pro- fessionally crafted, e ective communication, minor improve- ments possible in legend font size.” No further correction re- quired. Supplementary Note 3: Scenario 3 extended protocol This section documents the complete interaction sequence for Scenario 3 (ERA5/MOSAiC Lagrangian validation), including verbatim user prompts,data characteristics, the resource- ecient retrieval strategy, and all autonomous error-recovery episodes. 3.1 Data and scienti c context The user loaded a published minute-resolution underway dataset from the rst leg of the MO- SAiC expedition (PS122/1; Dataset DOI: 10.1594/PANGAEA.935221). The dataset contains 5 120,719 records across 22 columns, with temporal coverage spanning 2019-09-20 to 2019-12-13 and positions ranging from 69.68 ◦ N to 86.60 ◦ N and 19.00 ◦ E to 138.14 ◦ E. 3.2 User prompts Two user prompts drove the work ow: (1) validate ERA5 againstship measurements while minimizing download, and (2) generate a Wind Rose diagram with 16 cardinal bins and Beaufort color coding. 3.3 Self-correction episodes Two autonomous error-recovery events occurred: Episode 1: ERA5 bounding-box edge e ect.After initial interpolation, 83 of 2,013 co-located points returned NaN. The agent expanded the retrieval domain by∼0.3 ◦ and re- downloaded the subset, producing zero NaN values. Episode 2: Wind rose bin-label mismatch.The Visualization Agent de ned 13 Beau- fort speed-bin edges but provided 13 labels (12 required). The agent parsed theValueError traceback, removed the extraneous label, and re-executed successfully. Table S4: Annotated execution log for Scenario 3 (ERA5/MOSAiC Lagrangian validation). Bold text indicates autonomous agent de- cisions. The work ow comprises two user prompts (Steps 1, 11) and fourteen autonomous agent actions, including two error-recovery episodes (Steps 6–8, bounding-box edge e ect; Steps 13–14, bin- label mismatch) and two re exive visualization quality control cy- cles (Steps 10, 16). # ActorActionContent / Reasoning 1UserPrompt“Validate ERA5 reanalysis against the ship’s measurements while minimizing data download size. Calculate the exact spa- tial bounding box and time range from the ship’s data. Retrieve ERA52m temperatureonly for that speci c region. Interpolate the grid to the ship’s hourly positions, convert to Celsius. Cre- ate a time-series plot comparing in-situ air temperature (column T) in blue against ERA5 temperature in red.” 2SupervisorPlanningDecomposes into three sub-tasks: (i) spatiotemporal bounds ex- traction and hourly resampling→DataFrame Agent; (i) ERA5 cloud-native retrieval and Lagrangian interpolation→Oceanog- rapher Agent; (i) time-series visualization→Visualization Agent. 3DataFrameTemporal infer- ence Parses theDate/Timecolumn and computes spatiotemporal bounds: bounding box [18.99 ◦ E–138.14 ◦ E, 69.68 ◦ N–86.60 ◦ N], time range 2019-09-20 to 2019-12-13. Resamples to hourly reso- lution (last observation per hour), producingN= 2,013 hourly positions. Saveshourly shiptable.csv. 4Oceanographer API callSubmits ERA5 cloud-native retrieval request for2m- temperaturewith bounding box [18.99 ◦ E–138.14 ◦ E, 69.68 ◦ N– 86.60 ◦ N] and temporal range matching the cruise. Receives Zarr store: 2,013 time steps×68 latitude×477 longitude (∼261 MB). 6 # ActorActionContent / Reasoning 5Oceanographer Code executionInterpolates ERA52mtemperature eld toN= 2,013 hourly ship positions using trilinear interpolation (xarray.interp, method=linear). Converts from Kelvin to Celsius. 6OceanographerQuality check (fail) Detects 83 of 2,013 interpolated values returnedNaN. Diagnoses the cause: ship positions near the bounding-box edge fall outside the ERA5 grid domain (ERA5 grid starts at 69.75 ◦ N / 19.00 ◦ E, but ship positions extend to 69.68 ◦ N / 18.99 ◦ E). 7OceanographerSelf-correctionExpands the retrieval domain by∼0.3 ◦ in all directions (69.4 ◦ N– 86.9 ◦ N, 18.7 ◦ E–138.45 ◦ E). Re-submits the ERA5 API request. Receives expanded Zarr store: 70 latitude×479 longitude grid. 8Oceanographer Code executionRe-interpolates using the expanded grid. Veri es: 0 NaN values across all 2,013 co-located points. Saves merged tableship - hourly TTTvsERA5t2m.csv. 9Visualization Code executionGenerates dual-line time-series plot: ship in-situ T (blue) vs. ERA5 interpolated 2 m temperature (red). Uses seaborn-whitegridstyle, 300 DPI. 10 VisualizationRe ection loopRe exive quality control scores the time-series gure at9/10: “Axis labels clearly visible, good contrast, no overlapping ele- ments.” No further correction required. 11 UserPrompt“Generate a Wind Rose diagram to visualize the prevailing wind patterns for the whole cruise track. Split the wind directions into 16 cardinal bins and color-code the wedges based on the Beaufort Wind Force scale.” 12 Visualization Code executionLoads raw dataset. Identi es true wind columns:d(direction, ◦ ) andff(speed, m s −1 ). Filters ton= 120,616 valid observa- tions. 13 Visualization Code execution (fail) De nes 13 Beaufort speed-bin edges ([0, 0.5, 1.6, ..., 28.5,∞]) but provides 13 category labels.pd.cutraises ValueError: “Bin labels must be one fewer than the number of bin edges.” 14 VisualizationSelf-correctionParses the traceback. Recognizes the o -by-one mismatch: 13 edges require 12 labels. Removes the extraneous “Hurricane (12)” label (not represented in data; max speed = 23 m s −1 ). Re-executespd.cutsuccessfully. Constructs 16-sector×12- category frequency table. 15 Visualization Code executionGenerates stacked polar bar chart (wind rose) using South Po- lar Stereographic convention (0 ◦ = North, clockwise). Applies viridiscolormap across Beaufort categories. Exports PNG (400 DPI) and SVG. 16 VisualizationRe ection loopRe exive quality control scores the wind rose at8/10: “Clear, legible, standard meteorological convention adhered to.” No further correction required. Supplementary Note 4: Scenario 4 extended protocol This section documents the complete interaction sequence for Scenario 4 (Biodiversity and Water Mass Coupling), including verbatim user prompts, statistical outputs, and autonomous error-recovery episodes. 7 4.1 Data and scienti c context The user loaded a biological dataset of jelly sh abundance collected during the RVG.O. Sars cruise (May 2013) in the North Atlantic (DOI: 10.1594/PANGAEA.829702). The dataset comprises 94 sampling stations. 4.2 User prompts Three progressive prompts drove the analysis: (1) divide stations into Eastern/Western tran- sects at−20 ◦ longitude; (2) calculate Shannon-Wiener Index and generatea violin plot; (3) re- trieve Copernicus Marineθ/Sdata and create a T-S diagram colored byAglantha digitale abundance. 4.3 Self-correction episodes Two autonomous error-recovery events occurred: Episode 1: Mathematical safety.During Shannon Index calculation, the agent en- countered aRuntimeWarning: divide by zero. It refactored the code to apply a mask (where(p>0)). Episode 2: Re exive visual re nement.The initial violin plot scored 3/10. The agent iteratively refactored across ve cycles (3→4→4→5→9) until scienti c standards were met. 4.4 Statistical summary Welch’s t-test:t=−3.39,p= 0.0014. MeanH ′ West = 0.81 (n= 34, Irminger Sea),H ′ East = 0.35 (n= 60, Norwegian Sea). Table S5: Annotated execution log for Scenario 4 (G.O. Sars jelly- sh transect analysis with Copernicus validation). Bold text indi- cates autonomous agent decisions. The work ow comprises three user prompts (Steps 1, 6, 20) and twenty-three autonomous agent actions, including two error-recovery episodes (Steps 7–8, divide- by-zero in log; Steps 11–19, extended re exive visualization quality control with ve iterations, score progression: 3→4→4→5→9/10), and two single-pass re ection cycles (Steps 5, 26). # ActorActionContent / Reasoning Phase 1: Geographical transect division 1UserPrompt“Divide the sampling stations into two geographical transects based on longitude to separate the Norwegian Sea leg from the Irminger Sea leg. Create two new dataframes: ‘Eastern Tran- sect’ and ‘Western Transect’ for stations, use−20 ◦ longitude as a border. Visualize these transects at the end.” 2DataFrameData engineer- ing Identi esLongitudeandEventcolumns. Splits at−20 ◦ E: Eastern Transect(60 rows, 18 unique stations) andWest- ern Transect(34 rows, 11 unique stations). 3DataFrameCode executionGenerates lon/lat scatter plot with cruise track, dashed border at−20 ◦ , and “Irminger Sea leg” / “Norwegian Sea leg” anno- tations. Savestransects maplonlatborder-20.png. 8 # ActorActionContent / Reasoning 4Visualization Code executionParsesDate/Time, splits transects, sorts by time. Regenerates the transect map with connected cruise-track lines, styled scat- ter points, and text annotations. 5VisualizationRe ection loopRe exive quality control scores the transect map at8/10: “Axis labels clear, colors distinguishable, no overlapping elements.” Accepted without further correction. Phase 2: Shannon–Wiener diversity analysis 6UserPrompt“Calculate the Shannon–Wiener Diversity Index for everysta- tion using all available taxa columns. Then, create a split violin plot (or box plot with jitter points) to statistically compare the biodiversity levels between the ‘Eastern’ and ‘Western’ tran- sects.” 7DataFrameEcological com- putation Computes station-level Shannon–WienerH ′ =− ∑ p i lnp i across 64 taxa. EncountersRuntimeWarning: divide by zero in logfrom zero-abundance proportions. 8DataFrameSelf-correctionRefactors computation: converts proportions to NumPy array, appliesnp.where(p>0, p*np.log(p), 0.0)to avoid log(0). Veri es: minimumH ′ = 0 (no negative values). 9DataFrameStatistical infer- ence Computes inferential statistics: Welch’st=−3.39,p= 0.0014; Mann–WhitneyU= 638,p= 0.0019; Cohen’sd=−0.81; Hedges’g=−0.81; rank-biserialr= 0.37. Western transect: ̄ H ′ = 0.81±0.71 vs. Eastern: ̄ H ′ = 0.35±0.47. 10 Visualization Code executionGenerates violin + jitter-strip plot(seaborn) with quartile in- ner lines,n-annotations, and signi cance bracket with Welch p-value. 11 VisualizationRe ection loop (iter. 1) Quality control scores3/10: “Critical failure:p-value annota- tion andn-labels overlap with violin body and data points.” 12 VisualizationSelf-correctionIncreases Y-axis headroom (ylim×1.45), repositionsn-labels and bracket above violins, adds white text-box backgrounds. 13 VisualizationRe ection loop (iter. 2) Quality control scores4/10: “Annotationn=34 still overlaps with data distribution; Welch annotation overlaps with data plot area.” 14 VisualizationSelf-correctionSwitchestoaxes-coordinateannotations (transform=ax.transAxes) to guarantee non-overlap re- gardless of data range. Placesn-labels at (0.25, 0.98) and (0.75, 0.98),p-value at (0.5, 1.08). Setstitle pad=28. 15 VisualizationRe ection loop (iter. 3) Quality control scores4/10: “Title text overlaps with itself, impacting readability.” 16 VisualizationSelf-correctionIncreasestitle padfrom 28 to 40 to separate the two-line title from thep-value annotation. 17 VisualizationRe ection loop (iter. 4) Quality control scores5/10: “Y-axis is missing units; critical for scienti c accuracy.” 18 VisualizationSelf-correctionAppends unit annotation to Y-axis label: “Shannon diversity, H ′ (unitless; natural log).” 19 VisualizationRe ection loop (iter. 5) Quality control scores9/10: “No critical failures. Text and labels clear, appropriately sized. Conventions met.” Accepted. Phase 3: Copernicus T–S diagram with species overlay 20 UserPrompt“Retrieve Copernicus Marine data for Potential Temperature (‘thetao’) and Salinity (‘so’) matching the cruise coordinates and depths. Create a Temperature–Salinity (T–S) Diagram. Color-code the scatter points by the abundance ofAglantha dig- itale(using a log scale).” 9 # ActorActionContent / Reasoning 21 OceanographerData pro lingLoads cruise data. ParsesDate/Time(UTC). Computes depth midpoints fromDepth top/Depth bot(range: 12.5– 900 m). Extracts spatiotemporal bounds: 59.38–68.79 ◦ N, −46.19–5.07 ◦ E, 2013-05-03 to 2013-05-20. 22 Oceanographer API callSubmits Copernicus Marine retrieval forcmemsmodglophy- my 0.083degP1D-m(GLORYS12V1): variablesthetao,so; do- main [−50 ◦ –10 ◦ E, 58–70 ◦ N, 0–1000 m]; temporal range 2013- 05-01 to 2013-05-22. Receives NetCDF: 22 days×35 depths× 145 lats×721 lons (∼1.3 GB). 23 OceanographerCo-locationImplements 4D nearest-grid-point matching using xarray.sel(method=’nearest’)across time (daily oor), latitude, longitude, and depth. Achieves 94/94 (100%) valid matchups. Savescruise samplesmatchedcmemsthetao- so.csv. 24 DataFrameCode executionLoads data, extractsA. digitaleabundance column. Com- putes log 10 (x+ 1) transformation. Produces tidy table: 29 unique sampling events. Savesaglantha digitale- abundance tidy.csv. 25 Visualization Code executionLoads both intermediate CSVs. Inner-joins matched CMEMS data (n= 94) withA. digitaleabundance (n= 29) on Eventkey. Generates T–S scatter plot: salinity onx-axis, potential temperature ony-axis, viridis colormap scaled to log 10 (abundance + 1). Adds colorbar, sample-count annotation (n= 94), and grid. 26 VisualizationRe ection loopQuality control scores9/10: “High professional quality. Axis la- bels clear, color gradient e ective, conventions correct, no over- lapping elements.” Accepted without further correction. Supplementary Note 5: Search benchmarking methodol- ogy 5.1 Retrieval architectures under evaluation We evaluated three retrieval architectures: (1) Baseline (standard Elasticsearch BM25 keyword matching); (2) Simple LLM (GPT-5.2 zero-shot query translation); (3) Agentic Search (our ReAct-based iterative agent). 5.2 Query dataset construction We constructed a dataset of 100 natural language queries strati ed across ve complexity categories (20 queries each): •Category 1: Speci c Entity Search.Requests for a known, unique dataset. •Category 2: Broad Thematic Search.General exploration without strict spatiotem- poral lters. •Category 3: Spatiotemporal Slicing.Queries requiring translation of semantic loca- tions to coordinates. •Category 4: Parameter-Speci c.Requests targeting speci c chemical or physical variables. 10 •Category 5: Cross-Domain Constraints.Complex queries combining biological, physical, and geological lters. 5.3 Evaluation protocol and metrics For each query, the top-5 retrieved datasets from each architecture were evaluated program- matically using Gemini 3 Pro Preview (Google DeepMind) as anautomated judge, scoring on a scale of 0–10 across ve semantic metrics: 1.Precision / Relevance (M1):Overall relevance to the scienti c question. 2.Spatial / Temporal Accuracy (M2):Adherence to the requested geographic and temporal constraints. 3.Parameter Coverage (M3):Whether the datasets contain the speci c variables re- quested. 4.Access / Usability (M4):Suitability for immediate computational analysis. 5.Noise Reduction (M5):Ability to lter out false positives. Final scores per architecture were computed as the arithmetic mean across all 100 queries for each metric. ScenarioErrors Error TypesVQC Cycles Recovery 1: Weddell Sea1API depth parameter validation11/1 2:HAUS- GARTEN 1Datetime type mismatch21/1 3:MO- SAiC/ERA5 2Bounding-box edge e ect, bin-label mis- match 22/2 4: Jelly sh T-S1Divide-by-zero in log61/1 Total5115/5 (100%) Table S6: Summary of autonomous error recovery across all validation scenarios. “Errors” counts distinct runtime failures resolved via kernel traceback introspection. “VQC Cycles” counts the total number of re exive visual quality-control(VQC) cycles performed, including both correction iterations (score below acceptance threshold) and single-pass evaluations. All errors were resolved without human intervention. Supplementary Note 6: Qualitative retrieval analysis To complement the aggregate benchmark statistics reportedin the main text, we present a qual- itative comparison of representative queries that expose distinct failure modes in the baseline architectures. Table S7 compares the top-ranked result returned by each method for four queries selected to span di erent categories of semantic complexity: temporal reasoning, measurement geometry disambiguation, geographic constraint enforcement, and scienti c product-level intent recognition. 11 6.1 Temporal reasoning (Q065) Query:“Salinity and temperature pro les from the Weddell Sea during winter 2013.” The term “winter” in a Southern Hemisphere context maps to June–September, a constraint that requires hemispheric awareness absent from keyword indices. The Baseline returned Ice- Bird Winter altimetry campaigns from 2019 and 2023, matchingthe keyword “Winter” while violating both the year and parameter constraints. The Simple LLM retrieved oceanographic pro les but selected summer cruises (ANT-XXIX/2, PS82), failingto invert the seasonal cal- endar. The Agentic Search correctly resolved the hemispheric constraint and ranked the ANT- XXIX/6 (AWECS) cruise, conducted June–August 2013, as the top result (Agentic mean: 9.8/10; Baseline mean: 1.8/10). 6.2 Source disambiguation (Q059) Query:“Satellite altimetry data for studying sea level variability.” This query exposes what we term thevalidation trap: archival metadata frequently contains phrases such as “used to validate satellite altimetry” in records describing in-situ mooring deployments. Both the Baseline and Simple LLM returned these ground-truth observations, confusing validation context with the product itself. The Agentic Search correctly distinguished between the two categories, prioritizing Level-3/4 gridded sea-surface height anomaly products from CryoSat-2, Envisat, ERS-2, and TOPEX/Poseidon. It was the only method to retrieve the foundational WOCE TOPEX/Poseidon datasets (Agentic mean: 9.6/10; Baseline mean: 1.0/10). 6.3 Geographic constraint enforcement (Q080) Query:“Ground-penetrating radar (GPR) data from a polar glacier.” The Baseline and Simple LLM returned GPR datasets from Alpineglaciers in Austria and Italy, matching “Glacier” and “Radar” while ignoring the “polar” constraint. The Agentic Search strictly enforced the geographic restriction: every result in its top-30 list originated from Antarctica or Greenland. It also demonstrated domain-awaresemantic expansion, recognising that “Radio-Echo Sounding” (RES) is the standard airborne equivalent of ground-based GPR in polar glaciology, thereby broadening recall without sacri cing precision (Agentic mean: 9.8/10; Baseline mean: 2.8/10). 6.4 Product vs. raw data intent (Q096) Query:“Coastal erosion rates from historical aerial photographsand LiDAR.” The Simple LLM and Baseline returned raw aerial imagery subsets (Rettelbach et al. 2024, occupying 17–18 positions in each result list), treating the input methodology as the retrieval target rather than the requested scienti c product (erosion rates). The Agentic Search pri- oritised derived products: all top-6 results were datasetsquantifying shoreline retreat rates (Yukon, Laptev Sea, Baltic Sea, Portugal), correctly interpreting the user’s intent as the out- put variable rather than the input data source (Agentic mean:9.0/10; Baseline mean: 2.6/10). 12 Query in- tent Baseline (Elasticsearch)Simple LLMAgentic Search (ours) Q065:Wed- dell Sea, Win- ter 2013 (tem- poral) Failure:ReturnsIceBird Winter 2019/2023altime- try. Wrong year, wrong pa- rameter. Partial:Returns CTD pro les but from Summer cruises(ANT-XXIX/2, PS82). Success:ReturnsANT- XXIX/6 (AWECS), the spe- ci c winter cruise (Jun–Aug 2013). Q059:Satel- lite altimetry (source) Failure:Returns lake water quality and biological tran- sects. Validation trap:Returns in-situ mooring data “used to validate” satellites. Success:Returns Level- 3/4 gridded SSHA from CryoSat-2,Envisat, TOPEX/Poseidon. Q080:Polar glacier GPR (spatial) Geo-drift:Returns GPR from Alpine glaciers (Aus- tria/Italy).Matches “Glacier+Radar”but ignores “Polar”. Mixed:Returns Antarc- tic data but includes Alpine mass balance stakes. Success:Strictly polar. All 30 results from Antarctica or Greenland. Expands to RES (Radio-Echo Sounding). Q096: Coastal ero- sionrates (intent) Failure:Returns vege- tation maps and oceano- graphic bottle data. Intent miss:Returns raw aerial imagery subsets (18/30 results from one campaign). Success:Returns derived shoreline retreat rates (Yukon, Laptev, Baltic, Portugal). Table S7: Head-to-head comparison of the top-ranked search result for four representative complex queries, illustrating qualitative failure modes in baseline methods that are resolved by the agentic architecture. Each query tests a distinct category of semantic reasoning: tempo- ral constraint resolution, source vs. validation data disambiguation, geographic ltering, and product-level intent recognition. Supplementary Note 7: Dynamic contextualization pipeline The Oceanographer Agent contextualizes local observationswithin global environmental elds through a three-stage pipeline. First, it extracts the spatiotemporal bounding box of the loaded observational dataset (latitude/longitude extremes, temporal range, depth coverage). Second, it maps the user’s scienti c request to the appropriate external product: ERA5 for atmospheric variables, GLORYS12V1 (Copernicus Marine Service) for ocean state variables, or satellite-derived products for surface elds such as chlorophyll-aor sea-surface temperature. Third, it constructs a minimal retrieval request, specifying only the required variable, pressure level or depth layer, and spatiotemporal subset, and executes the download via the provider’s API. For high-volume atmospheric data, the system accesses the ERA5 reanalysis through the Earthmover Arraylake repository, a cloud-native chunked Zarr store that supports lazy loading viaxarrayandfsspec. This architecture reduces data transfer from terabytes (full global elds) to the megabyte-scale subsets actually required for the analysis. Once retrieved, the agent performs geometric alignment: for Eulerian comparisons, it co-locates observations to the nearest grid point using the Haversine distance; for Lagrangian platforms (e.g., ship tracks), it performs 4D interpolation along the trajectory in time, latitude, longitude, and depth. 13