Metodo numerico efficiente per la modellazione della dinamica del ghiaccio marino

Jinlun Zhang
Centro di Scienza Polare, Università di Washington, Seattle
W. D. Hibler III
Scuola di Ingegneria Thayer, Dartmouth College, Hanover, New Hampshire

Riassunto:
Presentiamo un metodo numerico altamente efficiente per risolvere i modelli di dinamica del ghiaccio marino che utilizzano reologie visco-plastiche. Questo metodo si avvale di un approccio semi-implicito per separare le equazioni del momento del ghiaccio nelle direzioni x e y, ottenendo così migliori proprietà di convergenza rispetto alle equazioni accoppiate. Oltre a velocizzare le soluzioni che utilizzano tecniche di rilassamento puntuale, abbiamo trovato che l’uso di una tecnica di sovrarilassamento lineare successivo, abbinata a un metodo di soluzione di sistemi matriciali tridiagonali, garantisce una convergenza particolarmente rapida. Questo procedimento può essere applicato anche alle equazioni della dinamica del ghiaccio in coordinate curvilinee ortogonali, fornite qui in forma esplicita per il caso particolare delle coordinate sferiche.

1. Introduzione

Una caratteristica saliente del ghiaccio marino nelle regioni polari è il suo movimento quasi incessante. Nei modelli dinamici del ghiaccio marino, questo movimento è solitamente descritto attraverso equazioni del momento che considerano la copertura di ghiaccio come un continuum bidimensionale. È noto da tempo che l’interazione tra i ghiacci è un processo fisico complesso e fortemente non lineare. Recentemente, diverse reologie plastiche non lineari sono state impiegate [Coon et al., 1974; Pritchard et al., 1977; Hibler, 1979; Flato e Hibler, 1992; Ip et al., 1991] per modellare queste interazioni non lineari nelle equazioni del momento del ghiaccio marino. Il metodo visco-plastico proposto da Hibler nel 1979 per simulare il flusso plastico è diventato molto popolare, poiché permette di rappresentare una gamma di leggi costitutive plastiche relativamente realistiche e complesse, caratterizzate da una forte non linearità. L’idea fondamentale del metodo “visco-plastico” consiste nell’approssimare la parte rigida o elastica di un continuum plastico con uno stato di creep molto lento.

L’inclusione delle forze di interazione non lineare del ghiaccio nelle equazioni del momento aumenta significativamente la difficoltà di risoluzione [ad es., Hibler, 1979]. Questo perché l’adozione di una reologia plastica, resistente alla compressione, impedisce di combinare le due equazioni del momento per le componenti di velocità u e v in una singola equazione complessa di velocità, come si fa di solito nei modelli di strato limite [ad es., McPhee, 1986; Thorndike e Colony, 1982]. Sebbene le equazioni accoppiate, come quelle proposte da Hibler [1979], possano essere risolte attraverso metodi espliciti di avanzamento temporale che sono quasi completamente ottimizzabili su computer vettoriali o multiprocessori [ad es., Ip et al., 1991], il passo temporale deve essere molto piccolo per garantire la stabilità delle soluzioni. Pertanto, il vantaggio dell’ottimizzazione completa viene annullato dal breve passo temporale. Di conseguenza, i metodi espliciti risultano impraticabili tranne che per applicazioni molto limitate.

Un metodo semi-implicito ad alta efficienza computazionale è stato proposto da Hibler nel 1979 per risolvere le equazioni del movimento del ghiaccio, includendo una reologia plastica con una curva di resa ellittica. Questo metodo unisce uno schema di passo temporale di tipo Euler modificato con una procedura di rilassamento successivo puntuale. In ogni passo temporale, si trattano i termini non lineari e si procede al passo successivo. Questo richiede un intenso lavoro di rilassamento puntuale, eseguito due volte per ogni passo temporale, per risolvere le equazioni accoppiate e lineari implicite per le componenti della velocità del ghiaccio, u e v. Poiché nel metodo semi-implicito le viscosità di massa e di taglio sono linearizzate, la soluzione ottenuta ad ogni passo temporale rappresenta solo un’approssimazione del flusso plastico. Tuttavia, ripetendo diverse volte la procedura semi-implicita in ciascun passo temporale fisico, è possibile ottenere una soluzione effettivamente plastica. Recentemente, Ip, nel 1993, ha esteso questa procedura di base per includere una varietà di curve di resa plastica più complesse.

Nonostante sia efficace a bassa risoluzione, questo metodo di rilassamento successivo puntuale risulta computazionalmente inefficiente per griglie grandi ad alta risoluzione. Una delle ragioni di questa inefficienza è che si deve ricorrere a un rilassamento limitato o a un sovrarilassamento a causa della complessità delle equazioni accoppiate u e v, come sottolineato da Ames nel 1977. È stato inoltre osservato che, con il rilassamento puntuale, il tasso di convergenza dipende fortemente dal massimo parametro di viscosità o di “strisciamento” presente nella griglia. Ridurne il valore può accelerare la convergenza, ma anche portare a una soluzione del flusso plastico meno realistica. Recentemente, Oberhuber nel 1990, e successivamente Holland e altri nel 1993, hanno applicato una variazione del sovrarilassamento successivo lineare per risolvere le equazioni del movimento del ghiaccio visco-plastico, adattate a un sistema di coordinate sferiche. In questo metodo, si utilizza una variazione del rilassamento per fila, eseguendo una scansione riga per riga nella griglia a differenze finite del modello, e risolvendo direttamente e simultaneamente la velocità del ghiaccio u e v in ogni riga con un efficiente solutore di matrici tridiagonali. Sebbene questo metodo rappresenti un miglioramento rispetto al rilassamento puntuale, continua a focalizzarsi sulla risoluzione delle equazioni completamente accoppiate.

In questo studio, presentiamo un nuovo approccio che migliora notevolmente l’efficienza computazionale, indipendentemente dal metodo di soluzione numerica adottato. La strategia principale è quella di separare le equazioni per le componenti di velocità u e v, in modo tale che le equazioni separate presentino migliori proprietà di convergenza e consentano l’impiego di metodi più efficienti per risolvere le equazioni implicite. Il flusso plastico viene quindi ottenuto attraverso una semplice procedura di avanzamento temporale. Questo metodo si basa su un trattamento semi-implicito delle equazioni del momento ad ogni livello dello schema di passi temporali di Euler modificato. Come dimostrato più avanti, le equazioni implicite risultanti dalla separazione possono essere risolte in maniera particolarmente efficiente utilizzando i solutori di matrici tridiagonali, seguendo l’approccio di Oberhuber del 1990.

Il documento descrive il metodo e mostra la sua efficacia in termini di efficienza computazionale e plausibilità delle soluzioni. Analizziamo anche il tasso di convergenza verso il flusso plastico e come questo tasso dipenda dal parametro massimo di viscosità di creep. È importante notare che abbiamo esteso questo metodo per includere diverse curve di resa plastica [per esempio, Ip, 1993], oltre alla curva di resa ellittica usata nella maggior parte degli esempi trattati.

2. Descrizione del Metodo

Secondo Hibler [1979], il movimento del ghiaccio marino è determinato da un equilibrio delle forze in gioco. La velocità del ghiaccio è espressa attraverso un vettore che include componenti direzionali lungo gli assi x e y. Fattori come la massa del ghiaccio per unità di superficie, il parametro di Coriolis, l’accelerazione gravitazionale, e l’altezza dinamica della superficie marina, influenzano questo movimento. A questi si aggiungono lo stress causato dall’aria e la resistenza dell’acqua, che sono non lineari, oltre alle forze di interazione interna del ghiaccio.

Le forze dovute allo stress dell’aria e dell’acqua vengono calcolate usando i coefficienti di attrito specifici per l’aria e l’acqua, considerando le densità di entrambi e gli angoli di deviazione del vento e delle correnti oceaniche geostrofiche.

La forza risultante dall’interazione interna dei ghiacci è rappresentata tramite un tensore di stress. In un sistema isotropico, questo tensore è legato alla velocità di deformazione e alla resistenza del ghiaccio attraverso una relazione visco-plastica non lineare. La velocità di deformazione del ghiaccio viene determinata osservando come variano le componenti di velocità relative alle direzioni spaziali. La resistenza del ghiaccio, invece, dipende dalla sua compattezza e dallo spessore. I parametri relativi alla “viscosità”, che includono la viscosità volumetrica e quella di taglio, sono funzioni degli invarianti del tasso di deformazione del ghiaccio e della sua resistenza, raggiungendo un valore massimo di “strisciamento” quando il tasso di deformazione è particolarmente basso.L’obiettivo di questo approccio è quello di simulare un comportamento rigido attraverso uno stato di strisciamento molto lento. Per la curva di resa ellittica impiegata, le viscosità non lineari si distinguono tra loro per un fattore costante. Nelle simulazioni discusse in questo studio, abbiamo assunto che questo fattore e la resistenza del ghiaccio siano indipendenti dal tasso di deformazione. È possibile esplorare reologie plastiche più elaborate consentendo variazioni più complesse delle viscosità e facendo dipendere la resistenza del ghiaccio dagli invarianti di deformazione, come mostrato in studi precedenti. Attraverso una scelta accurata della dipendenza funzionale delle viscosità e della resistenza del ghiaccio dal tasso di deformazione, il metodo può essere esteso anche per adattarsi a curve di resa basate sulla frattura.

2.1. Schema Numerico

Per semplicità, il metodo numerico per risolvere le equazioni è inizialmente descritto utilizzando un sistema di coordinate rettangolari. Questo metodo, tuttavia, può essere applicato anche alle equazioni del momento in qualsiasi sistema di coordinate curvilinee ortogonali. Queste equazioni curvilinee verranno illustrate brevemente nella sezione successiva per evidenziare la loro similitudine strutturale.

Le componenti della forza di interazione del ghiaccio, derivanti dalla legge costitutiva, sono calcolate all’interno del sistema di coordinate cartesiane. Le equazioni che descrivono queste componenti possono essere organizzate e semplificate facilmente. Le formule per le velocità del ghiaccio e le viscosità sono funzioni degli invarianti del tasso di deformazione. È rilevante sottolineare che i termini di advezione non lineare sono stati esclusi dalla trattazione, in quanto molto meno significativi rispetto agli altri termini.

Nelle equazioni per le componenti di velocità u e v, non mostrate qui, tutti i termini sul lato destro sono calcolati utilizzando la velocità del ghiaccio già determinata, mentre le componenti della velocità sul lato sinistro sono aggiornate. Questo approccio disaccoppia efficacemente le equazioni per le componenti u e v, trasformando ogni equazione in una formula lineare implicita indipendente, risolvibile separatamente. Applicando il metodo delle differenze finite a ciascuna di esse, si ottiene un sistema di equazioni lineari simultanee, configurando una matrice simmetrica, definita positiva e con dominanza diagonale. Queste caratteristiche permettono l’uso di tecniche di sovrarilassamento senza il rischio di divergenza, risultando così in una rapida convergenza.

In particolare, il trattamento di una determinata equazione implica un’elaborazione implicita del termine di Coriolis, tenendo fissi i termini di interazione del ghiaccio, e può essere sostituito da un’altra equazione più semplice da risolvere attraverso calcoli algebrici. Tale equazione emerge dalla sottrazione di due specifiche equazioni. È importante notare che, in assenza di interazioni nel ghiaccio, questa procedura si semplifica in una soluzione completamente implicita delle equazioni di deriva libera. Nel caso in cui si voglia mantenere un smorzamento neutrale dei termini inerziali, il lato destro può essere modificato opportunamente.

Come dimostrato in un’apposita appendice, questa formulazione è altresì numericamente stabile.

Le Equazioni (3) e (4) possono essere risolte attraverso tecniche di rilassamento sia puntuale che lineare, impiegando un risolutore di matrici tridiagonali. Per un’efficienza ottimale, quando si adotta il rilassamento lineare, è cruciale applicare un sovrarilassamento successivo riga per riga (nella direzione x) per le equazioni componenti di u, e colonna per colonna (nella direzione y) per quelle di v. Questo perché il primo termine nelle equazioni della componente u enfatizza la direzione x, mentre il primo termine delle equazioni di v predilige la direzione y. È stato osservato che conducendo il sovrarilassamento per righe nelle equazioni della componente v, il tasso di convergenza può essere triplicato rispetto al rilassamento colonna per colonna.

L’Equazione (3) può essere ridefinita secondo un’equazione alle differenze finite che segue lo schema di Hibler del 1979, per una soluzione efficiente. Questo schema può essere implementato tramite una scansione sistema per sistema in ogni riga, utilizzando tecniche come l’algoritmo di Thomas. Similmente, le equazioni della componente v possono essere risolte con un approccio per colonne.

Dettagli importanti di questo algoritmo e delle sue variazioni per diverse condizioni al contorno sono disponibili nell’Appendice B.

Per un confronto, nella soluzione originaria di Hibler del 1979, le derivate spaziali e i termini di Coriolis erano trattati in modo implicito. Questo approccio è stato mantenuto anche nella metodologia di Oberhuber del 1990. È interessante notare che nel lavoro originale di Hibler, il calcolo del coefficiente di dragaggio non lineare dell’acqua non era temporizzato adeguatamente. Hibler e Ackley nel 1983 hanno però osservato che, in presenza di forze di interazione del ghiaccio minime, il trattamento non corretto di questi termini può generare soluzioni errate, a meno che non si adotti un approccio di centratura modificato, come quello proposto da Euler.

Si nota inoltre che, in situazioni specifiche dove la viscosità di taglio è assente, come nella soluzione visco-plastica dei fluidi cavitanti, il processo di scansione sequenziale delle equazioni non necessita di un approccio iterativo LSOR, ma basta una singola soluzione tridiagonale per risolvere le equazioni implicite. Questo rende il metodo particolarmente efficiente dal punto di vista computazionale per questo tipo di reologia del ghiaccio.

La soluzione suddivisa in tre passaggi non garantisce sempre una risposta plastica alla velocità del ghiaccio, che può essere necessaria in alcuni contesti. Per raggiungere una soluzione plastica, è possibile adottare dei “passi di pseudo tempo”, durante i quali le viscosità vengono aggiornate ad ogni intervallo. Con “passo di pseudo tempo” si intende una procedura in cui, dopo aver completato un intero passo temporale con la tecnica di Euler modificato, la velocità finale viene trattata come punto di partenza per tutti i termini, ad eccezione di quello inerziale, e il processo viene ripetuto mantenendo invariata la resistenza del ghiaccio. Un approccio ancora più efficiente prevede la ripetizione del solo secondo livello del passo temporale modificato di Euler. In questa fase, la velocità del ghiaccio, centralizzata, continua a essere utilizzata per calcolare la resistenza all’avanzamento dovuta all’acqua, ma è la velocità appena ottenuta (dal primo livello o dal pseudo tempo precedente) a essere impiegata per aggiornare le viscosità e, se necessario, anche la velocità centralizzata. Dopo vari cicli di pseudo tempo, la velocità del ghiaccio si avvicinerà a una soluzione plastica. È necessario un ulteriore passaggio di correzione implicita, realizzabile come descritto in precedenza, usando la velocità derivante dall’ultimo passo di pseudo tempo.

2.2. Formulazione in sistemi di coordinate curvilinee ortogonali

Questo metodo numerico si rivela efficace anche nella risoluzione delle equazioni del momento in sistemi di coordinate curvilinee ortogonali. La procedura base per adattare il modello a un sistema di coordinate curvilinee qualsiasi prevede di esprimere la forza del ghiaccio e la legge costitutiva in forma diadica, adatta a qualsiasi sistema di coordinate curvilinee ortogonali. Il tensore di stress viene definito in termini della velocità del ghiaccio. In un sistema curvilineo, le operazioni sui gradienti devono tenere conto delle unità vettoriali, che variano in funzione delle coordinate curvilinee. Qui, per esempio, queste operazioni sono analizzate specificamente per un sistema di coordinate sferiche. Considerando che il modello è bidimensionale, le componenti della forza di interazione del ghiaccio vengono calcolate per coordinate geostrofiche bidimensionali. Le coordinate geostrofiche bidimensionali si riferiscono a un sistema di coordinate sferiche sulla superficie terrestre, senza variazioni nella direzione radiale.

I risultati, basati sui lavori di Malvern del 1969 e sull’Appendice B, sono stati formulati considerando variabili come il raggio della Terra, la longitudine e la latitudine. Si considerano componenti specifiche in cui i coefficienti sono influenzati dagli invarianti del tensore di tasso di deformazione, rappresentati in forma diadica.

Utilizzando le coordinate sferiche, specifiche componenti sono calcolate prendendo in considerazione variabili cruciali come la longitudine e la latitudine. Applicando questo metodo, è possibile derivare e disaccoppiare le equazioni del momento in modo analogo a quanto fatto con le equazioni cartesiane. Questo permette l’uso dello stesso schema di avanzamento temporale a tre livelli già descritto per risolvere le equazioni del momento. Anche lo schema di differenze finite spaziali proposto da Hibler nel 1979 è applicabile a questo contesto.

I primi test effettuati utilizzando coordinate sferiche hanno dimostrato un incremento dell’efficienza numerica grazie al disaccoppiamento, simile a quello osservato utilizzando coordinate rettangolari.

3. Risultati delle Simulazioni

Per confrontare il metodo di disaccoppiamento con le precedenti tecniche di soluzione, abbiamo condotto una serie di simulazioni del ghiaccio usando un computer parallelo vettoriale Alliant FX-80. Abbiamo analizzato metodi di rilassamento puntuale per entrambe le formulazioni delle equazioni. Per le equazioni disaccoppiate, abbiamo anche testato un metodo di sovrarilassamento successivo combinato con un risolutore di matrici tridiagonali. Nei metodi per le equazioni disaccoppiate, questi sono identificati come DPSOR (rilassamento puntuale disaccoppiato) e DLSOR (sovrarilassamento successivo disaccoppiato), mentre per le equazioni accoppiate abbiamo utilizzato PSOR (rilassamento puntuale accoppiato). Abbiamo utilizzato la griglia di Arakawa B, posizionando le velocità agli angoli delle celle della griglia e le viscosità al centro, in tutti i casi. Questa griglia agevola l’integrazione con i modelli di circolazione oceanica.

Questi tre metodi sono stati applicati a due modelli di ghiaccio con diverse risoluzioni che coprono l’Artico, il Mare di Groenlandia e il Mare di Norvegia, come illustrato nella Tabella 1. A parte la risoluzione, i modelli di ghiaccio sono simili al modello di Hibler del 1979 e includono un collegamento con equazioni esplicite per la conservazione della concentrazione e dello spessore del ghiaccio. La Figura 1 mostra la configurazione della griglia del modello di ghiaccio con risoluzione di 40 km. Per pilotare le simulazioni del ghiaccio, sono stati usati i dati giornalieri dei campi di forzatura atmosferica del 1981. Dettagli più approfonditi sui campi di forzatura sono stati forniti da Hibler e Zhang nel 1993.

la Tabella 1 dettaglia due modelli di ghiaccio utilizzati nelle simulazioni:

1. Modello 1:
   • Ha una risoluzione di griglia di 80 chilometri e copre una griglia di 65 per 66 celle.
   • Il modello utilizza un intervallo di tempo di un giorno per ciascun passo di calcolo.
   • Le simulazioni con questo modello durano per un periodo totale di 30 giorni.
2. Modello 2:
   • Questo modello ha una risoluzione più fine con 40 chilometri per cella e una griglia più grande di 130 per 102 celle.
   • Il passo temporale usato è di mezza giornata, permettendo simulazioni più dettagliate nel tempo.
   • La durata totale delle simulazioni per questo modello è di 10 giorni.

Inoltre, la tabella descrive come la resistenza del ghiaccio nei modelli sia determinata in base a vari fattori come la massa media di ghiaccio per unità di superficie e la sua densità o concentrazione. Questi fattori sono cruciali per capire come il ghiaccio reagirà a diverse condizioni ambientali e come verrà modellata la sua dinamica nelle simulazioni.

La Figura 2 mostra un grafico che confronta la media dell’energia totale rispetto agli intervalli di tempo di simulazione nel modello di ghiaccio con una risoluzione di 40 km. Il grafico presenta tre diverse linee, ciascuna rappresentante un metodo differente di simulazione:

1. Linea tratteggiata: Questa linea mostra l’energia totale media calcolata usando il metodo PSOR, che è un tipo di procedura di rilassamento.
2. Linea continua: Rappresenta l’energia totale media calcolata con un metodo standard, non specificato nel dettaglio.
3. Linea punteggiata: Indica la differenza nelle velocità calcolate tra il metodo standard e un nuovo metodo.

Sull’asse orizzontale del grafico sono rappresentati gli intervalli di tempo dei passaggi della simulazione, che vanno da 2 giorni a 1/8 di giorno. Si può notare che riducendo l’intervallo di tempo tra i passaggi, l’energia totale media tende a diminuire sia per il metodo PSOR che per il metodo standard. Anche la differenza nelle velocità tra i metodi calcolati mostra una diminuzione.

Il grafico è utile per vedere come intervalli di tempo più brevi nelle simulazioni possano migliorare la precisione dei risultati, riducendo sia l’energia totale, che può indicare potenziali errori o instabilità, sia le discrepanze nelle velocità calcolate dai diversi metodi. Questo aiuta a capire l’efficacia dei diversi metodi di simulazione sotto condizioni di passi temporali più frequenti.

3.1. Confronti dei Risultati Numerici

Salvo diversa indicazione, abbiamo confrontato il metodo di rilassamento puntuale per le equazioni accoppiate (PSOR) con il metodo di rilassamento lineare per le equazioni disaccoppiate (DLSOR). Tuttavia, un altro metodo chiamato DPSOR ha prodotto risultati quasi identici per le equazioni disaccoppiate. Per tutti i metodi è stato utilizzato un comune parametro di sovrarilassamento impostato a 1.5, e la loro sensibilità a questo parametro sarà discussa in seguito.

Per valutare il comportamento generale del campo di velocità del ghiaccio simulato, la Figura 2 illustra l’energia totale media e la differenza media totale di velocità tra i metodi, calcolate su un periodo di integrazione di un mese in funzione della lunghezza dei passi temporali per il modello di ghiaccio di 40 km. L’energia totale è calcolata come la somma dei quadrati delle velocità del ghiaccio in ciascuna cella della griglia, mentre la differenza totale di velocità è la somma dei quadrati delle differenze di velocità tra i due metodi. Questo dimostra che, con la riduzione degli intervalli dei passi temporali, i trattamenti semi-impliciti diversi tendono a convergere verso la stessa soluzione. Si osserva che i valori di energia totale ottenuti dai due metodi si avvicinano e la discrepanza nella velocità si riduce con la diminuzione dell’intervallo dei passi temporali. Passi temporali di mezza giornata sono stati sufficienti per mantenere le discrepanze entro un intervallo contenuto.

Entrambi i metodi tentano di approssimare il flusso plastico tramite passaggi temporali impliciti, quindi non è possibile preferirne uno in modo assoluto. Di conseguenza, un modo più incisivo per valutare l’efficacia dei due metodi potrebbe essere il confronto del loro comportamento quando si avvicinano a soluzioni di equilibrio per la velocità del ghiaccio plastico. A questo scopo, è stato realizzato uno studio di sensibilità utilizzando condizioni di vento e correnti oceaniche costanti nel tempo ma variabili nello spazio. Per analizzare la soluzione di equilibrio plastico, sono stati mantenuti costanti anche lo spessore e la compattezza del ghiaccio. La velocità iniziale del ghiaccio è stata impostata a zero, e i venti sono stati attivati al primo passo temporale.

La Figura 3 mostra le variazioni totali della velocità rispetto allo stato di equilibrio e l’energia totale, rappresentate in funzione del tempo. Si osserva che, una volta attivati i campi di forzatura, l’energia totale calcolata con entrambi i metodi converge rapidamente verso la soluzione di equilibrio plastica. È stato anche confermato che il metodo DLSOR raggiunge il flusso plastico più rapidamente rispetto ad altre curve di resa plastica, come dimostrato da Ip nel 1993. Un altro aspetto significativo mostrato dalla Figura 3 è che con l’aumentare della tolleranza, i risultati ottenuti tramite il metodo di rilassamento puntuale tendono a raggiungere l’equilibrio plastico più lentamente. La tolleranza qui è definita come la massima differenza ammessa tra i valori vecchi e nuovi della velocità del ghiaccio al termine di un ciclo iterativo. Di conseguenza, una tolleranza più bassa non garantisce una soluzione precisa delle equazioni linearizzate. Tuttavia, utilizzando il rilassamento lineare, la soluzione tridiagonale integra le condizioni al contorno in una soluzione esatta di parte delle equazioni differenziali, riducendo così la sensibilità alla tolleranza in questo scenario.

Infine, è importante notare che lo studio di sensibilità discusso può essere esteso al caso speciale di una reologia di fluido cavitante, come descritto da Flato e Hibler nel 1992. Mantenendo inalterati tutti gli altri parametri, la viscosità di taglio è stata impostata a zero, trasformando così la reologia visco-plastica di Hibler del 1979 in un fluido cavitante, ovvero un fluido privo di resistenza al taglio. Questa formulazione del fluido cavitante differisce leggermente da quella proposta da Flato e Hibler, poiché la viscosità di volume e di conseguenza lo stress compressivo sono leggermente influenzati dalla deformazione di taglio, come indicato da Ip nel 1993. Tuttavia, queste comparazioni rimangono utili per valutare l’efficacia dei metodi numerici. La Figura 4 illustra come il nuovo metodo raggiunga più rapidamente la soluzione di equilibrio, evidenziando che il disaccoppiamento delle equazioni nel nuovo approccio permette una soluzione esatta delle equazioni separate, utilizzando metodi di soluzione tridiagonale senza la necessità di ulteriori rilassamenti.

La Figura 1 illustra la configurazione della griglia per un modello di ghiaccio con risoluzione di 40 km, usato per studiare le dinamiche del ghiaccio nelle regioni dell’Artico, della Groenlandia e dei mari Norvegesi. La griglia copre vasti territori che includono parti del nord della Russia, l’Alaska, e la Groenlandia, estendendosi fino ai mari Norvegesi.

Ogni punto sulla griglia rappresenta un’area di terra o mare di 40 km per 40 km, permettendo un’analisi dettagliata e accurata del comportamento del ghiaccio marino. La densità di questa griglia assicura una risoluzione spaziale adeguata per modellare processi complessi come le interazioni tra il ghiaccio marino e le correnti oceaniche, oltre agli effetti delle variazioni climatiche su queste aree sensibili.

Questa configurazione di griglia è fondamentale per le simulazioni che mirano a comprendere e prevedere i cambiamenti nel clima e l’impatto di tali cambiamenti sulle regioni polari, fornendo uno strumento essenziale per i ricercatori nel campo della climatologia e degli studi ambientali.

La Figura 3 mostra due grafici principali per ogni uno dei tre pannelli etichettati come a, b e c, ognuno dei quali esplora i risultati ottenuti con diversi livelli di tolleranza nella convergenza dei calcoli.

Nella parte superiore di ogni pannello, è rappresentata l’energia totale del modello di ghiaccio. L’energia totale riflette l’attività dinamica del ghiaccio, ovvero quanto sia energico e in movimento in un determinato momento.

Nella parte inferiore di ogni pannello, è mostrata la deviazione della velocità del ghiaccio rispetto a uno stato ideale di equilibrio, calcolata su tutti i punti della griglia che contengono ghiaccio marino. Questa misura aiuta a capire quanto si discostano le velocità attuali del ghiaccio da uno stato ideale e stabile.

Le linee continue nei grafici rappresentano i risultati ottenuti usando il metodo PSOR, un tipo di calcolo che aggiorna la posizione del ghiaccio punto per punto. Le linee tratteggiate mostrano i risultati del metodo DLSOR, un approccio che aggiorna le informazioni lungo linee intere della griglia, potenzialmente offrendo un modo più rapido per raggiungere l’equilibrio.

I tre pannelli si distinguono per i livelli di tolleranza utilizzati nei calcoli:

• Pannello a: Impiega la tolleranza più stretta, mirando a una maggiore precisione.
• Pannello b: Utilizza una tolleranza standard, bilanciando precisione e velocità di calcolo.
• Pannello c: Adotta una tolleranza più ampia, potenzialmente accelerando i calcoli ma con una precisione ridotta.

Questi grafici sono utili per valutare come i differenti metodi e le varie tolleranze influenzano la rapidità e l’accuratezza con cui i modelli di ghiaccio raggiungono uno stato di equilibrio, offrendo intuizioni preziose su quale configurazione potrebbe essere più efficace in contesti di ricerca o operativi reali.

3.2. Confronti dell’Efficienza Computazionale

Entrambi i metodi simulano in modo simile il flusso plastico, quindi l’attenzione si focalizza principalmente sull’efficienza computazionale del metodo di disaccoppiamento. Per valutarla, abbiamo analizzato il numero di iterazioni richieste da ciascun metodo e il tempo di calcolo consumato per risolvere la velocità del ghiaccio, utilizzando come base un elemento computazionale del computer Alliant FX80, senza ricorrere a tecniche di parallelizzazione.

Un aspetto cruciale di tutti questi metodi di soluzione è quanto possono essere ottimizzati attraverso una scelta adeguata del parametro di sovrarilassamento. Abbiamo esplorato questa questione conducendo una serie di test per ottenere soluzioni plastiche di equilibrio con il modello da 40 km, variando il parametro di sovrarilassamento ω. Questi test, che prevedevano una simulazione di 10 giorni partendo da un campo di velocità iniziale nullo, sebbene poco efficienti, hanno permesso di comparare i diversi tassi di convergenza dei metodi.

La Figura 5a mostra l’impatto di vari parametri di sovrarilassamento sull’efficienza computazionale, rappresentando il numero medio di iterazioni in funzione del parametro di sovrarilassamento per ciascun metodo. Grazie alle caratteristiche superiori di convergenza delle equazioni disaccoppiate, è stato possibile utilizzare un parametro di sovrarilassamento notevolmente superiore rispetto a quello utilizzabile con le equazioni accoppiate, che non sono stabili per valori di sovrarilassamento oltre 1.5. I confronti dell’efficacia relativa di ciascun metodo, mostrati nella Figura 5b, rivelano che il metodo di sovrarilassamento lineare successivo disaccoppiato ha un punto di efficienza massima intorno a ω = 1.7, rendendolo circa cinque volte più efficiente del miglior risultato ottenuto con il metodo di rilassamento puntuale disaccoppiato, il quale ha un punto di efficienza massima più stretto intorno a ω = 1.9.

Per un test più pratico dell’efficienza computazionale dei vari metodi, abbiamo utilizzato il campo di velocità del passo temporale precedente come valore iniziale. Questo approccio è comune nei procedimenti di avanzamento temporale e aiuta a valutare l’efficacia operativa dei diversi metodi.

Per valutare questo metodo in confronto all’approccio che parte da una stima iniziale nulla, abbiamo realizzato simulazioni sui modelli di ghiaccio con risoluzioni di 80 e 40 km, impiegando un parametro di sovrarilassamento di 1.5 per tutti i modelli. I risultati sono illustrati nelle Figure 6 e 7. Come emerge dalle figure, l’utilizzo di una buona stima iniziale (Figura 7) fornisce risultati simili a quelli ottenuti con l’approccio di partenza da zero (Figura 6), ma con un numero significativamente inferiore di iterazioni e, di conseguenza, un minor consumo di tempo di calcolo.

Il metodo di disaccoppiamento si dimostra ancora estremamente più efficiente sia in termini di numero di iterazioni che di tempo di calcolo, anche se ora il metodo di rilassamento puntuale disaccoppiato risulta essere solo marginalmente più efficiente rispetto a quello accoppiato. Tuttavia, secondo quanto mostrato dalla Figura 5, questa situazione cambierebbe se venisse utilizzato un valore di sovrarilassamento più elevato. Il rilassamento lineare disaccoppiato rimane comunque il metodo di soluzione più efficace.

Abbiamo anche cercato di valutare in modo limitato gli effetti della vettorizzazione e della concorrenzializzazione sull’efficienza computazionale dei vari metodi. Questo è stato fatto confrontando le simulazioni che utilizzavano un numero variabile di elementi computazionali su un computer parallelo vettoriale Alliant FX80. I risultati sono presentati nella Figura 8. In ogni caso, il grado di ottimizzazione riflette la capacità del compilatore di vettorializzare e gestire in parallelo i codici, il che può non corrispondere all’efficienza massima ottenibile.

La Figura 4 mostra un grafico che rappresenta la deviazione totale della velocità rispetto allo stato di equilibrio in un modello di ghiaccio con risoluzione di 80 km, che utilizza una particolare reologia denominata “fluido cavitante”. Questo grafico mette a confronto due metodi di rilassamento per vedere quanto velocemente ciascuno di essi riesce a stabilizzare la velocità del ghiaccio verso un equilibrio desiderato.

• Linea continua (PSOR): Questa linea, che rappresenta il metodo di rilassamento puntuale (Point Successive Overrelaxation), mostra una diminuzione graduale della deviazione della velocità. Questo indica che il metodo sta lavorando per portare la velocità del ghiaccio verso lo stato di equilibrio, ma con una velocità di convergenza moderata.
• Linea tratteggiata (DLSOR): La linea tratteggiata rappresenta il metodo DLSOR (Decoupled Line Successive Overrelaxation), un approccio di rilassamento lineare disaccoppiato. Questo metodo mostra una riduzione più rapida e decisa della deviazione rispetto al PSOR, suggerendo una convergenza più veloce verso lo stato di equilibrio.

Il grafico dimostra che entrambi i metodi sono efficaci nel ridurre la deviazione della velocità, ma il metodo DLSOR è visibilmente più efficiente, raggiungendo una condizione di equilibrio più rapidamente rispetto al metodo PSOR. Questo è particolarmente rilevante in scenari di modellazione dove la rapidità di stabilizzazione può avere implicazioni importanti sulla precisione e sull’utilità delle simulazioni condotte.

La Figura 5 presenta due grafici che analizzano l’efficienza dei vari metodi numerici usati nel modello di ghiaccio con risoluzione di 40 km.

Grafico (a): Questo pannello mostra il numero medio di iterazioni richieste in funzione del parametro di sovrarilassamento per tre differenti procedure numeriche:

• Quadrati: Indicano il metodo di rilassamento puntuale per le equazioni accoppiate (PSOR).
• Diamanti: Rappresentano il metodo di rilassamento lineare per le equazioni disaccoppiate (DLSOR).
• Triangoli: Denotano il metodo di rilassamento puntuale per le equazioni disaccoppiate (DPSOR).

Il grafico rivela che l’aumento del parametro di sovrarilassamento tende a ridurre il numero di iterazioni necessarie, in particolare per il metodo di rilassamento lineare (DLSOR), che mostra una diminuzione più pronunciata rispetto agli altri.

Grafico (b): Questo pannello visualizza il miglioramento relativo frazionario di ogni metodo rispetto a un approccio di rilassamento diretto. Si misura la percentuale di iterazioni risparmiate rispetto a un approccio più semplice e diretto, indicando quanto ciascun metodo sia più efficiente.

• I quadrati, diamanti e triangoli mostrano rispettivamente l’efficacia del PSOR, del DLSOR e del DPSOR.

Questo grafico dimostra che il metodo DLSOR (diamanti) è generalmente il più efficiente, specialmente a valori più elevati del parametro di sovrarilassamento, suggerendo che questo metodo può ottenere risultati comparabili con meno iterazioni, migliorando notevolmente l’efficienza computazionale.

Insieme, questi grafici offrono una visione approfondita di come i vari metodi di rilassamento impattino il tempo e l’efficacia computazionale delle simulazioni del ghiaccio marino, fornendo indicazioni utili su quale metodo potrebbe essere più vantaggioso per specifiche configurazioni di simulazione.

La Figura 6 confronta l’efficienza computazionale di tre metodi di rilassamento iterativo applicati a due differenti modelli di simulazione del ghiaccio, misurata sia in termini di tempo di CPU (secondi) che di numero medio di iterazioni necessarie per completare le simulazioni.

Analisi per modello:

1. Modello 1:
   • PSOR (rilassamento puntuale accoppiato): Mostra il tempo di CPU più elevato e il maggior numero di iterazioni, indicando un’efficienza relativamente bassa.
   • DPSOR (rilassamento puntuale disaccoppiato): Riduce sia il tempo di CPU che il numero di iterazioni rispetto al PSOR, mostrando un miglioramento nell’efficienza.
   • DLSOR (rilassamento lineare disaccoppiato): Risulta essere il metodo più efficiente per il Modello 1, minimizzando ulteriormente il tempo di CPU e il numero di iterazioni.
2. Modello 2:
   • PSOR: Analogamente al Modello 1, questo metodo richiede il maggiore impegno computazionale sia in termini di tempo che di iterazioni.
   • DPSOR: Presenta un’efficienza migliore rispetto al PSOR ma non è il più ottimale.
   • DLSOR: Dimostra nuovamente di essere il metodo più efficiente, con il minor tempo di CPU e il minor numero di iterazioni necessarie.

Osservazioni generali: I metodi disaccoppiati, DPSOR e DLSOR, mostrano un’efficacia superiore in entrambi i modelli rispetto al tradizionale PSOR. In particolare, il metodo DLSOR emerge come il più performante, offrendo significative riduzioni nel tempo di CPU e nel numero di iterazioni richieste. Questo suggerisce che la tecnica di rilassamento lineare disaccoppiato può essere estremamente vantaggiosa per migliorare l’efficienza computazionale nelle simulazioni del ghiaccio.

Questi risultati sottolineano l’importanza della scelta del metodo di rilassamento nella progettazione di simulazioni computazionali efficienti, specialmente in contesti dove il tempo di calcolo e le risorse sono fattori critici.

La Figura 7 fornisce un confronto dell’efficienza computazionale tra tre metodi di rilassamento iterativo, applicati a due diversi modelli di simulazione del ghiaccio, il Modello 1 e il Modello 2. Questa efficienza è valutata sia in termini di tempo di CPU consumato che del numero medio di iterazioni necessarie per ciascun metodo.

Analisi dettagliata:

• Barre bianche: Indicano il tempo di CPU impiegato da ciascun metodo, misurato in secondi.
• Barre a scacchi: Rappresentano il numero medio di iterazioni richieste da ciascun metodo per completare la simulazione.

Metodi Valutati:

• PSOR (Point Successive Overrelaxation): Metodo di rilassamento puntuale accoppiato.
• DPSOR (Decoupled Point Successive Overrelaxation): Metodo di rilassamento puntuale disaccoppiato.
• DLSOR (Decoupled Line Successive Overrelaxation): Metodo di rilassamento lineare disaccoppiato.

Risultati:

1. Modello 1:
   • Il PSOR richiede il maggior tempo di CPU e il maggior numero di iterazioni.
   • Il DPSOR mostra un miglioramento significativo rispetto al PSOR, con una riduzione del tempo di CPU e del numero di iterazioni.
   • Il DLSOR emerge come il metodo più efficiente, minimizzando sia il tempo di CPU che il numero di iterazioni.
2. Modello 2:
   • Simili tendenze si osservano nel Modello 2, con il PSOR che registra ancora il tempo di CPU più alto e il maggior numero di iterazioni.
   • Il DPSOR e il DLSOR migliorano entrambi l’efficienza, con il DLSOR che si conferma come il più efficace, riducendo drasticamente sia il tempo di CPU che il numero di iterazioni rispetto agli altri metodi.

Considerazioni aggiuntive: In questi test, la stima iniziale per la soluzione di ogni procedura iterativa utilizza il campo di velocità ottenuto al passo temporale precedente, permettendo potenzialmente una convergenza più rapida verso la soluzione finale e un uso più efficiente delle risorse computazionali.

In conclusione, la Figura 7 illustra chiaramente i vantaggi dei metodi di rilassamento disaccoppiati, in particolare del DLSOR, che supera gli altri metodi in termini di efficienza computazionale, sia per il tempo di calcolo che per il numero di iterazioni necessarie. Questo fornisce un insight prezioso sull’ottimizzazione dei processi di simulazione in contesti computazionali intensivi.

Anche se altre tecniche di soluzione tridiagonale potrebbero essere ottimizzate in modo più efficiente, i risultati attuali offrono comunque diversi spunti importanti. In particolare, la vettorizzazione e la concorrenzializzazione migliorano la velocità delle soluzioni solo del 30% per i metodi di rilassamento puntuale, indicando che la maggior parte del tempo di calcolo è consumato dalla parte della procedura che non può essere vettorializzata. Il miglioramento relativo dovuto alla vettorizzazione è maggiore nel metodo DLSOR, suggerendo che questo metodo contiene meno elementi non vettorializzabili. Inoltre, il metodo DLSOR beneficia anche di una maggiore velocità grazie all’uso di processori paralleli, probabilmente perché le equazioni u e v vengono risolte simultaneamente su differenti unità computazionali. Teoricamente, un miglioramento simile potrebbe essere possibile anche per il metodo DPSOR con un computer più ottimizzato.

Si osserva inoltre che, con le equazioni disaccoppiate, il numero di operazioni richieste durante ogni ciclo di rilassamento è leggermente inferiore rispetto a quello necessario con le equazioni accoppiate. Un’efficace strategia per incrementare l’efficienza dei metodi di rilassamento puntuale consiste nel ridurre la viscosità massima di “creep”. Tuttavia, questo può portare a soluzioni molto diverse dal flusso plastico desiderato, un risultato indesiderabile. Questo accade perché nel metodo visco-plastico lo stato di creep viscoso è utilizzato per approssimare il flusso plastico rigido. Sono stati condotti esperimenti numerici per valutare l’impatto della viscosità massima sulle prestazioni di questi metodi in termini di efficienza computazionale. In generale, la grandezza di questa viscosità determina l’efficacia nella simulazione del flusso plastico rigido, con valori più elevati che ne migliorano l’approssimazione e valori più bassi che degradano la reologia del ghiaccio a un comportamento viscoso lineare. Il valore standard utilizzato nei lavori di Hibler è sufficientemente alto per replicare adeguatamente il flusso plastico e ha dimostrato di essere efficace nelle analisi di lungo periodo in cui il ghiaccio rimane sostanzialmente immobile.L’impatto del valore massimo di viscosità sui due metodi è rappresentato nella Figura 9. Qui, il tempo di CPU utilizzato per la dinamica del ghiaccio in un’integrazione mensile del modello ghiaccio-oceano con una risoluzione di 80 km (effettuata su un singolo elemento computazionale) è graficato rispetto al logaritmo del valore massimo di viscosità diviso il suo valore standard. Dal grafico si osserva che il metodo PSOR è fortemente influenzato da questo valore di viscosità, mentre il metodo DLSOR mostra una dipendenza molto più debole. Questa caratteristica si presume sia dovuta alla tecnica di rilassamento lineare successivo utilizzata dal DLSOR. Tuttavia, il vantaggio del metodo DLSOR in termini di efficienza computazionale rispetto al metodo PSOR tende a diminuire progressivamente con la degradazione della reologia del ghiaccio verso un comportamento viscoso lineare.

La Figura 8 presenta un confronto del tempo di CPU, misurato in secondi, per tre diversi metodi di rilassamento iterativo impiegati nel “Model 1”. Questo confronto è strutturato per mostrare l’efficacia dell’ottimizzazione computazionale attraverso la vettorizzazione e l’uso di più elementi computazionali (CEs).

Configurazioni di calcolo:

• Barre piene: Tempo di CPU su un singolo elemento computazionale senza vettorizzazione.
• Barre a scacchi: Tempo di CPU su un singolo elemento computazionale con ottimizzazione completa.
• Barre vuote: Tempo di CPU quando il calcolo è distribuito su tre elementi computazionali con ottimizzazione completa.

Analisi dei metodi:

• PSOR (Point Successive Overrelaxation): Registra i tempi di CPU più alti in tutte le configurazioni, indicando che è il metodo meno efficiente tra quelli testati.
• DPSOR (Decoupled Point Successive Overrelaxation): Mostra un miglioramento nell’efficienza rispetto al PSOR, con riduzioni nel tempo di CPU in tutte le configurazioni.
• DLSOR (Decoupled Line Successive Overrelaxation): Dimostra di essere il metodo più efficiente, soprattutto quando ottimizzato e utilizzato su tre elementi computazionali, evidenziando una significativa diminuzione del tempo di CPU.

Osservazioni generali:

• L’impiego della vettorizzazione e di più elementi computazionali porta a una notevole riduzione del tempo di CPU per tutti i metodi, ma l’effetto è particolarmente marcato per il DLSOR.
• L’efficienza aumenta notevolmente con l’uso di tre CEs ottimizzati, sottolineando l’importanza di sfruttare le capacità di concorrenzializzazione e vettorizzazione per migliorare l’efficienza computazionale nelle simulazioni di modelli di ghiaccio.

Questi risultati mettono in luce come la scelta dell’hardware e delle tecniche di ottimizzazione possano influenzare significativamente le prestazioni dei metodi computazionali, e come tali tecniche possano essere strategicamente utilizzate per ottimizzare i processi computazionali in scenari di simulazione intensiva.

La Tabella 2 presenta il tempo di CPU, espresso in secondi, necessario per calcolare una soluzione plastica utilizzando il metodo DLSOR (Decoupled Line Successive Overrelaxation). La simulazione è stata condotta su un modello di ghiaccio con un passo temporale giornaliero per un periodo di un mese.

Il metodo impiegato, noto come “pseudo time stepping”, è utilizzato per migliorare la convergenza del modello durante l’integrazione. La tabella specifica il tempo di CPU richiesto per diversi numeri di “pseudo time steps”:

• A 1 pseudo time step, il tempo di CPU è il più basso, pari a 2027 secondi.
• Incrementando i pseudo time steps a 5, il tempo di CPU sale a 4287 secondi.
• Con ulteriori aumenti, il tempo di CPU continua a crescere, raggiungendo 6018 secondi a 10 pseudo time steps e 9449 secondi a 30 pseudo time steps.

L’aumento del tempo di CPU con il numero di pseudo time steps indica che un maggiore numero di aggiornamenti della soluzione, effettuati durante l’integrazione mensile, richiede più risorse computazionali. Questo comporta un trade-off tra la precisione e la rapidità della simulazione: più passaggi intermedi possono potenzialmente migliorare la qualità della simulazione ma a spese di tempi di calcolo più lunghi.

La Figura 9 visualizza il tempo di CPU necessario per la soluzione della dinamica del ghiaccio in un modello di risoluzione di 80 km, analizzato per un’integrazione mensile. Il grafico mette in relazione questo tempo di CPU con il logaritmo del rapporto tra la viscosità massima utilizzata nel modello e un valore di viscosità standard. Questo confronto è stato effettuato utilizzando tre diversi metodi numerici: PSOR, DPSOR e DLSOR, rappresentati rispettivamente da una linea continua, una tratteggiata e una tratteggiata-puntata.

Dettagli chiave del grafico:

• Asse orizzontale (X): Mostra il logaritmo del rapporto di viscosità, indicando come variazioni nella viscosità influenzino il processo computazionale.
• Asse verticale (Y): Rappresenta il tempo di CPU in secondi, fornendo una misura diretta dell’efficienza computazionale.

Osservazioni dai dati:

• Metodo PSOR: Questa linea mostra che il tempo di CPU aumenta notevolmente con l’aumento del rapporto di viscosità, indicando una forte sensibilità di questo metodo alle variazioni di viscosità.
• Metodo DPSOR: Con un incremento meno marcato del tempo di CPU rispetto al PSOR, questa linea suggerisce che il DPSOR è moderatamente più efficiente e meno sensibile alle variazioni di viscosità.
• Metodo DLSOR: Questa linea dimostra la minore sensibilità ai cambiamenti nel parametro di viscosità, evidenziando la sua stabilità e efficienza superiore rispetto agli altri due metodi.

In conclusione, la Figura 9 evidenzia l’efficienza comparativa dei tre metodi di simulazione rispetto alle variazioni di un parametro chiave come la viscosità. Il metodo DLSOR emerge come il più efficiente e il meno influenzato dalle variazioni della viscosità, rendendolo una scelta potenzialmente migliore per simulazioni in condizioni variabili.

3.3. Efficienza Computazionale del Pseudo Time Stepping

Il Pseudo Time Stepping si rivela un meccanismo efficiente per raggiungere un flusso completamente plastico a ogni passo temporale. La sua efficacia si deve principalmente al fatto che non necessita del completo utilizzo del passo temporale di Euler modificato e che, avvicinandosi al flusso plastico, il numero di iterazioni richieste per ogni pseudo time step diminuisce. Per dimostrare l’efficienza di questo metodo, abbiamo analizzato il consumo di CPU per ottenere una soluzione plastica utilizzando il Pseudo Time Stepping nel modello 1, applicando una procedura di soluzione tridiagonale disaccoppiata. La Tabella 2 illustra il tempo di CPU impiegato con vari numeri di pseudo time steps. È importante notare che il consumo di CPU non aumenta in modo proporzionale all’aumentare dei pseudo time steps; ciò è dovuto al fatto che le modifiche alla soluzione diventano meno frequenti con l’incremento dei pseudo time steps.

In pratica, utilizzando circa 15 pseudo time steps, le soluzioni ottenute si avvicinano notevolmente al flusso plastico. Questo è evidenziato nella Figura 10, dove abbiamo rappresentato gli stati di stress normalizzati per alcuni punti casuali della griglia computazionale del modello 2 per 1, 5, 15 e 100 pseudo time steps. Come si può osservare, già cinque pseudo time steps offrono una buona approssimazione al flusso plastico per i passi temporali fisici di 1/4 di giorno usati in questo studio. Questo risultato è stato confermato anche dall’analisi dell’energia totale in funzione del tempo. Infatti, nella pratica è quasi impossibile distinguere tra le serie temporali di energia ottenute con cinque e cento pseudo time steps.

La Figura 10 illustra come variano gli stati di stress in un modello di ghiaccio, utilizzando diversi numeri di pseudo time steps per ogni scenario di simulazione con passi temporali di un quarto di giorno. I grafici rappresentano quattro differenti configurazioni: con 1, 5, 15 e 100 pseudo time steps.

Dettagli dei Grafici:

• In ogni grafico, i punti neri rappresentano la misurazione degli stati di stress in punti specifici della griglia computazionale, scelti casualmente. Questi punti sono stati raccolti dal passo temporale numero 40, e uno ogni dodici punti è stato selezionato per la visualizzazione.

Osservazioni Chiave:

• Nel grafico con 1 pseudo time step, i punti mostrano una grande dispersione, indicando una varietà ampia nelle misurazioni di stress.
• Aumentando i pseudo time steps nei grafici successivi (5 e 15), si osserva una progressiva concentrazione dei punti attorno a una curva più definita, suggerendo una stabilizzazione degli stati di stress nel modello.
• Nel grafico con 100 pseudo time steps, i punti si allineano strettamente lungo una traiettoria chiara, dimostrando che la soluzione si è stabilizzata notevolmente e riflette un comportamento più coeso e prevedibile degli stati di stress.

Conclusione: Questi grafici dimostrano efficacemente che aumentando il numero di pseudo time steps, il modello raggiunge una rappresentazione più stabile e accurata del comportamento del ghiaccio, con una variabilità ridotta negli stati di stress osservati. Questo suggerisce che per ottenere simulazioni più precise del flusso plastico nel ghiaccio, l’utilizzo di un maggior numero di pseudo time steps può essere particolarmente vantaggioso.

4. Osservazioni Conclusive

Un nuovo metodo numerico per la risoluzione dei modelli di ghiaccio marino visco-plastico è stato presentato e si è dimostrato molto più efficiente rispetto ai precedenti metodi di rilassamento puntuale. Questo approccio è ampiamente applicabile e, sebbene sia stato utilizzato qui per una specifica curva di resa plastica, si adatta efficacemente a diverse curve di resa plastica, come evidenziato da Ip nel 1993. È stato provato che il metodo produce soluzioni stabili che concordano bene con le soluzioni ottenute risolvendo le equazioni accoppiate con metodi di rilassamento puntuale, come riportato da Hibler nel 1979. Questa accuratezza si mantiene anche utilizzando intervalli di tempo non eccessivamente ridotti. Le soluzioni tipiche derivanti da questo metodo evidenziano importanti vantaggi: si avvicinano rapidamente a una vera soluzione di equilibrio plastico e sono relativamente poco sensibili alla precisione richiesta per la soluzione di rilassamento.

Il fulcro di questo metodo è il disaccoppiamento semi-implicito delle equazioni del momento del ghiaccio. Questo trattamento separa le equazioni del momento del ghiaccio marino in componenti u e v, migliorando le proprietà di convergenza delle equazioni implicite rimanenti e permettendo l’uso di metodi iterativi più efficaci. L’impiego di un risolutore tridiagonale, combinato con una tecnica di rilassamento lineare, ha mostrato una convergenza particolarmente rapida. Un altro punto di forza nell’utilizzo dei risolutori tridiagonali è che l’aumento del consumo di CPU è meno marcato rispetto ai metodi di rilassamento puntuale, anche quando la viscosità massima cresce.

La combinazione di queste caratteristiche ha portato a una significativa riduzione dei tempi di calcolo, rendendo i modelli di ghiaccio marino visco-plastico praticamente più gestibili per griglie ad alta risoluzione di grande dimensione o per simulazioni climatologiche globali. La procedura di time-stepping, che aggiorna le viscosità non lineari, offre anche una base conveniente per modellare flussi plastici completi mediante l’uso di vari pseudo-time steps. In pratica, questo metodo permette di modellare flussi plastici completi con solo un incremento da tre a quattro volte del tempo di calcolo. Di conseguenza, questo metodo apre la strada all’accesso a un’ampia varietà di reologie del ghiaccio plastico oltre a quelle esaminate in questo studio.

Appendice A: Analisi della Stabilità

Per illustrare le caratteristiche fondamentali di stabilità del processo di time-stepping, è utile esaminare un caso particolare delle equazioni del momento in cui sia le viscosità di taglio che quelle volumetriche sono costanti nello spazio. Inoltre, per semplicità, si considerano costanti anche i coefficienti di resistenza idrodinamica dell’acqua. Con queste semplificazioni, le equazioni del momento si presentano in una forma che evidenzia come il legame principale tra le equazioni u e v nell’interazione del ghiaccio sia rappresentato dalla viscosità volumetrica.

L’approccio principale consiste nel trasferire alcuni dei termini viscosi sul lato destro dell’equazione, dove vengono trattati in modo esplicito. La questione chiave è quindi determinare se restano sufficienti termini impliciti per assicurare la stabilità del sistema.

Il metodo base di time-stepping si avvale di un passo temporale di tipo Euler modificato per affrontare i termini non lineari. Successivamente, viene eseguito un passo correttivo implicito per risolvere in modo implicito i termini relativi alla forza di Coriolis e alla resistenza idrodinamica che non sono sulla diagonale principale. Questo ultimo passo temporale implicito rende le equazioni completamente implicite nel caso non ci sia interazione tra i ghiacci.

In assenza di termini non lineari nelle equazioni (A1) e (A2), il processo necessita solo di un passo temporale diretto seguito da un passo correttivo. Procediamo quindi con un’analisi di stabilità usando il metodo di Fourier, considerando una componente specifica di Fourier per le variabili u e v, e determinando se un eventuale fattore di amplificazione supera l’unità. Dato il nostro interesse per la stabilità con passi temporali molto grandi, assumiamo che l’intervallo di tempo tenda all’infinito, così da poter annullare il termine derivativo del tempo.

Adottando il metodo di Fourier, ipotizziamo forme specifiche per u e v basate su esponenti complessi. Sotto queste ipotesi, trasformiamo le equazioni (A1) e (A2) per isolare certi coefficienti. La soluzione di queste equazioni segue un approccio a due fasi: inizia con un passo provvisorio seguito da un passo correttivo implicito. Durante questo passo correttivo, combiniamo certi coefficienti in una relazione specifica dove la somma di due coefficienti particolari risulta essere uno.

Questo processo dimostra come, anche con un’applicazione rigorosa del time-stepping, il sistema mantenga caratteristiche di stabilità essenziali, garantendo l’efficacia del metodo anche con l’incremento della complessità del modello.

Concetto di base nel passaggio implicito

Il passo implicito di base consiste nel mantenere i valori dei termini viscosi utilizzati inizialmente, mentre si risolvono implicitamente tutti i rimanenti termini di resistenza idrodinamica e le forze di Coriolis. Per dimostrare quando è necessario il passo correttivo implicito per le forze di Coriolis e i termini di resistenza idrodinamica fuori diagonale, concentriamoci su due equazioni specifiche. Queste equazioni producono una matrice che amplifica la velocità del vettore (u, v), con autovalori definiti dalle soluzioni di queste equazioni.

Se gli autovalori risultano avere magnitudini minori di uno, la soluzione sarà stabile. Tuttavia, se la grandezza di certi coefficienti supera quella di altri e i valori di certe variabili sono abbastanza piccoli, ciò può condurre a una situazione di instabilità.

Nell’analisi di stabilità generale, possiamo eliminare i termini u e v dalle equazioni. Risolvendole per il vettore velocità (u, v), otteniamo espressioni aggiornate per il passo temporale successivo.

Per rendere ancora più chiare le espressioni degli autovalori della matrice di amplificazione, consideriamo il caso in cui un particolare coefficiente sia piccolo, che rappresenta il caso più critico. In questo scenario, i numeratori delle espressioni nelle equazioni si semplificano, producendo nuove relazioni per u e v nel passo temporale seguente.

Nelle equazioni indicate, abbiamo preservato i termini Cd nel numeratore laddove necessario per assicurare che il numeratore resti finito quando i valori di kx e ky si avvicinano a zero. È evidente che i fattori di amplificazione per le variabili u e v risultano essere inferiori a uno quando il parametro 0 è minore di 1/2. Tuttavia, quando kx e ky si riducono ulteriormente, i fattori di amplificazione restano inferiori a uno a condizione che 0 sia minore di 1/2 e il valore di Cd sia abbastanza piccolo.

Per confermare il caso in cui il valore di Cj è piccolo, possiamo ricorrere a una soluzione esatta degli autovalori della matrice di amplificazione. Impostando le magnitudini degli autovalori uguale a uno e analizzando l’equazione, troviamo che i valori possono raggiungere o superare la magnitudine uno solo se Cd è inferiore a un certo valore, con 0 minore di 1/2, e questa condizione si verifica quando kx e ky sono estremamente bassi, il che rappresenta la situazione più instabile.

Questo processo aiuta a comprendere l’importanza di gestire accuratamente i termini viscosi e di resistenza idrodinamica nelle simulazioni numeriche, evidenziando come e quando i vari termini e coefficienti influenzino la stabilità delle soluzioni del modello.

Appendice B: Rilassamento Lineare con l’Algoritmo di Thomas

Nel contesto di una specifica equazione per determinare la componente della velocità x, u, lungo una riga, intendiamo risolvere per la velocità x nei vari punti della riga. Questa equazione collega la velocità x in ciascuna colonna a un termine noto che può dipendere dalle velocità u in altre righe o dalle componenti di velocità v.

Con le condizioni al contorno di Dirichlet stabilite, miriamo a risolvere questa serie di equazioni lineari per i valori che vanno dal secondo all’ultimo, adattandoci alle condizioni al contorno che specificano le velocità agli estremi della riga.

Per semplificare, possiamo assumere che la velocità all’estremità della riga sia zero, modificando di conseguenza il termine noto. In questo caso, i coefficienti dell’equazione sono definiti per gli indici da 2 a n – 1. Per risolvere questo gruppo di equazioni, impieghiamo l’Algoritmo di Thomas, un metodo iterativo di eliminazione gaussiana che necessita di due passaggi lungo la riga: un passaggio in avanti e uno in ritorno. Questa tecnica permette di trattare efficacemente e sistematicamente le equazioni lineari lungo una singola riga, garantendo una soluzione accurata e coerente con le condizioni al contorno impostate.

Durante il passaggio in avanti del nostro processo, stabiliamo due costanti indicizzate, g_i e f_i, iniziando con le prime e aggiornando successivamente le altre per ogni indice dal terzo fino all’ultimo meno uno. Questo metodo ci permette di calcolare fi utilizzando i valori precedentemente determinati.

Per ottenere le soluzioni per la velocità u, applichiamo un metodo che procede all’indietro, partendo dall’ultimo elemento fino al secondo, utilizzando le condizioni al contorno stabilite.

Quando implementiamo il rilassamento per righe, sostituiamo i valori originali di u con nuovi valori calcolati mediante un parametro di sovrarilassamento, scelto per ottimizzare l’efficienza. Questi nuovi valori sono derivati dalle fasi di calcolo precedenti. Procediamo poi a risolvere la riga successiva in modo analogo.

Questo approccio può essere adattato anche per le componenti della velocità y, utilizzando però il rilassamento per colonne. In circostanze dove il parametro di sovrarilassamento è unitario, è possibile applicare anche un rilassamento per blocchi.

La condizione al contorno standard prevede una velocità zero ai confini terrestri, una scelta che, grazie alla plasticità non lineare del modello, facilita un naturale scorrimento discontinuo in questi punti. I bordi liberi del ghiaccio marino sono gestiti dinamicamente impostando la resistenza del ghiaccio a zero vicino ai confini terrestri artificiali. Inoltre, un comportamento di scorrimento libero può essere implementato vicino ai confini terrestri utilizzando la formulazione di fluido cavitante discussa nel testo.

https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/96JC03744