Ajatusmallini metsien hakkuiden vaikutuksesta ilmastomuutokseen

Webissä käyty keskustelu metsien vaikutuksesta ilmastonmuutokseen herätti ajatuksen: "Tässäpä oiva ilmiö pähkäiltäväksi!" Ilmiön tutkimisen voisi aloittaa lukemalla kaiken olennaisen metsien kasvusta ja ilmastonmuutoksesta. Vaan mistä tietäisin, mikä on minulle olennaista? Pitäisikö lukea kaikki aiheista kirjoitettu? Minusta tulisi ehkä maailman ainoa, joka tietäisi kaiken metsistä ja ilmastonmuutoksesta, vaan riittäisikö ikäni moiseen luku-urakkaan. Osaisinko, jaksaisinko tehdä mitään kaikella tiedollani? Kun kaikkitietävänä alkaisin jakaa tietoani, kiinnostaisiko ketään enää? Ehkä ajat ja ongelmat olisivat vallan toiset.

Lienee parempi ensin miettiä, mitä ilmiöstä tietää ja sen jälkeen hakea tilkettä tietojensa aukkoihin sen mukaan kuin tarvetta ilmenee. Voisi jutella kaverien kanssa, kysyä opettajalta. Joukossa tieto jalostuu, jos porukalla on rakentavan vuoropuhelun taito. Harmi kyllä minulle on tarjolla juttukaveriksi vain vaimoni, jolla on omat kiireensä.

**Jos tuntuu, ettei asiasta mitään tiedä, kannattaa lähteä kyselemään opastusta. Osaamistaan voi levittää tietojensa ja taitojensa lähialueille, mutta on ajan hukkaa yrittää ponnistaa liikkeelle tyhjän päältä. Tosin moni vain luulee ettei tiedä mitään, kun ei ole varma, onko oma tieto oikeaa virallista tietoa. Tottahan kaikki jotain tietävät; puut kasvavat, hakkuissa puut tapetaan ja viedään metsästä, aurinko lämmittää maapalloa jne. Kiperämpi kysymys saattaa olla, miten metsien kasvu ja ilmaston lämpeneminen liittyvät toisiinsa. Kasvihuoneilmiöistä puhutaan, mutta eihän metsät kasvihuoneessa kasva. Pitää kysyä kaverilta, opettajalta tai kokeiltava, mitä hakukone kertoo kasvihuoneilmiöstä. Hakukoneen löydökset saattavat olla ymmärrettäviä tai pitää kysyä apua, millaisesta hausta olisi apua.

Käsityksensä ilmiöstä kannattaa laittaa muotoon, jonka muutkin ymmärtävät niin on selkeämpi jutella aiheesta. Ymmärrys ilmiöstä saattaa jalostua ja täydentyä.

Metsien kasvu ja ilmastonmuutos ovat dynaamisia ilmiöitä täynnä keskinäisvaikutuksia. Niinpä päätin laittaa oman käsitykseni siitä differentiaaliyhtälöiksi, jotka ratkaisen tietokoneella. Tällaisten simulaatioiden tulokset ovat yleensä havainnollisia ja tuovat ilmiöistä esille piirteitä, joita ei muuten välttämättä huomaisi.

Differentiaaliyhtälöiden numeerinen ratkaisu on tärkeä osa monia videopelejä, ilmastonmuutoksen tutkimusta, lentokoneiden suunnittelua yms. Harmi, ettei sitä taideta koulussa edes mainita.

Koulussa harjoitustehtävät on yleensä valmiiksi räätälöity tietyn yksittäisen asian harjoitteluun niin, että ratkaisuun on yksi, suora valmis polku, josta ei ole syytä poiketa. Elävän ilmiöiden tarkastelu rönsyilee väkisinkin, eikä ole niinkään selvää, mitkä polut vievät kelvolliseen ratkaisuun tai mikä yleensä on kelvollinen ratkaisu. Tämäkin juttu rönsyilee muistutukseksi siitä, että rönsyily on olennainen osa ilmiön tutkimista. Elävään elämään kuuluu myös pähkäily, onko tässä mitään järkeä, koulussa riittää tarkistaa, että vastaus on oikea.

Webin keskustelupalstoilla tulee esille monenlaisia näkökohtia, mutta kommenttien perusteella on vaikea hahmottaa, miten eri keskustelijat näkevät kokonaisuuden ja miten heidän kommenttinsa siihen liittyy. On myös vaikea hahmottaa, mistä loppujen lopuksi ollaan eri mieltä, mistä samaa mieltä ja mitkä olisivat olennaiset tarkemmin selvitettävät seikat.

Keskusteltaessa metsien hakkuiden vaikutuksesta ilmastonmuutokseen on esitetty mm. seuraavanlaisia näkemyksiä:

• Hiilinieluja tarvitaan myöhemminkin, joten hakkuiden kasvattaminen muutamaksi vuosikymmeneksi ei haittaa;
• Suomen metsät ovat niin nuoria ja kasvun alussa, että ilmastonmuutoksen kannalta järkevintä olisi antaa niiden järeytyä mieluummin hillitsemällä hakkuita kuin niitä lisäämällä;
• Metsiä kannattaa hakata, että kasvuvauhti pysyy hyvänä;
• 'Kypsä' metsä ei toimi hiilinieluna;
• 'Kypsä' metsä on hiilivarasto;
• Koska puusta tehty biodiesel on uusiutuvaa, se on hiilineutraali polttoaine, joten siihen kannattaa investoida;
• Puunkäytön ilmastovaikutukset riippuvat siitä, mitä hakatulle puulle tehdään;
• Hoidettu metsä sopeutuu paremmin ilmastonmuutokseen;
• Riittää, että metsä kasvaa enemmän kuin sitä hakataan;
• Hakatun metsän tilalle istutetaan uutta (joten hakkuilla ei ole vaikutusta?);
• Hiilidioksidi edistää metsän kasvua;
• Suomen ilmaston lämpeneminen edistää metsiemme kasvua;
• Polttoaineeksi käytetty hake on kasvaessaan jo sitonut poltossa vapautuvan hiilidioksidin;

Useimmat väitteista ovat tosia, melkein tosia tai ainakin puolitotuuksia, mutta monen osalta jää hämäräksi sen kuviteltu merkitys ilmaston lämpenemiselle. Osaa väitteistä käytetään tavalla, joka pistää epäilemään, että vaikutusmekanismeista on erilaisia ajatusmalleja. Tuntuu, että ilmiöiden ajallista käyttäytymistä on vaikea hahmottaa. Tuntuu, että ilmiötä kokonaisuutena ei ole mietitty loppuun saakka.

Fast thinking

Usein webissä keskustellaan Kahnemanin 'fast-thinking'-moodissa: Esitetään väitteitä oman intuition — omien ajatusmallien — pohjalta ja mahdollisesti poimitaan jostain todisteeksi omaan intuitioon sopiva lehtijuttu. Usein keskustelun käynnistääkin jonkun omaa intuitiota tukeva klikkiotsikko. Selvityksistä saatetaan poimia joku luku oman väitteen perusteluksi selvittämättä, miten tulokseen on päädytty, mistä se oikeasti kertoo, kuinka tarkasti ja mikä luvun merkitys on kokonaisuuden kannalta.

Näin käy meille kaikille. Ei auta, vaikka tietää näin olevan.

Keskustelu jää pinnalliseksi, ellemme ennakkoluulottomasti halua tutustua ajatusmalleihin ( = käyttöteorioihin), joihin eri keskustelijat perustavat väitteensä.

Omia ajatusmallejaan voi yrittää havainnollista vertauksin. Voi esimerkiksi verrata ylimääräistä hiilidioksidia ilmakehässä veteen veneen pohjalla? "Uusiutuvan biodieselin käyttäminen on kuin äyskäröisi vettä keulasta perään. OK, jos vettä on vain vähän pohjalla, mutta ei riitä, jos vene on uppoamaisillaan. Toki parempi kuin fossiilisten polttoaineiden käyttäminen, mikä vastaa veden äyskäröimistä järvestä veneeseen."

Vertauksilla voi onnistua vaikuttamaan mielipiteisiin, mutta vertaukset voivat johtaa harhaan, eivät todista mitään ja voivat johtaa keskustelun aivan sivuraiteille.

Slow-thinking

Keskustelua voisi mielestäni yrittää syventää havainnollistamalla omaa ajatteluaan kaaviokuvien ja kuvaajien avulla. Toisten olisi helpompi perustella näkemyksiään minulle, kun ymmärtäisivät ajatteluani ja näkisivät sen puutteet.

Koska ihmisen aivot eivät ole mikään numeronmurskaaja, päättelymme ja selityksemme fysikaalisista ilmiöistä perustuu kvalitatiivisiin malleihin. Sanotaan, että nuori metsä kasvaa vanhaa paremmin, hoidettu metsä sitoo hyvin hiiltä toisin kuin hoitamaton jne. Kvalitatiivisessa päättelyssä löydetään raja-arvoja — 'ääripäitä': "Hiljan istutettu taimikko ei juurikaan sido hiiltä", "Vanha metsä ei enää sido lisää hiiltä". Yksi tärkeä raja-arvo tuntuu olevan nolla: Laitetaan painoa sille, ollaanko hiilen sidonnassa plussalla vai miinuksella. Ainakin itse olen kysellyt, millä strategialla hiiltä saataisiin eniten varastoitua — joko metsiin, rakennuksiin tms.

Kuvitelmieni tarkastelua

Seuraavassa havainnollistan yksinkertaisten simulaatioiden avulla käsitystäni metsän kasvusta ja hakkuiden vaikutuksesta ilmaston lämpenemiseen. Kuvailen siis ajatusmalliani — yksinkertaistettua kuvitelmaani todellisuudesta. Tämä juttu ei siis todista mitään todellisuudesta.

Sovellan hyvin yksinkertaista matemaattista mallinnusta. Käyttämäni menetelmä ei millään tavalla muistuta ilmastomallinnuksessa käytettyä lähestymistapaa

Tarkastelu saattaa kuitenkin osoittaa ajatusmalleissani sisäisiä ristiriitoja ja ristiriitoja todellisuuden kanssa. Joku muu nähnee siinä helposti virheitä ja puutteita tai toteaa minun keskittyneen toisarvoisiin seikkoihin.

Webistä olisin ehkä löytänyt ohjelmiston, jolla olisin voinut paremmin tarkastella metsien hakkuun vaikutusta ilmastoon. Tarkoitukseni on kuitenkin esitellä python-ohjelmoinnin soveltamista yksinkertaisiin matemaattisiin tehtäviin, joten koodasin kaiken itse.

Selitän ohjelmoinnin periaatteita niin, että tästä voi olla hupia ohjelmoinnista kiinnostuneille, vaikka metsien hakkuut eivät kiinnostaisikaan. (Selittely tietysti tympii ohjelmointia tuntevia.)

Blogikirjoituksessa hakkuun vaikutus metsan hiilensidontaan on mallintamiani seikkoja selitetty havainnollisesti. Eihän tuokaan ottaen todista mitään, mutta auttanee keskustelun jäsentämisessä.

Ohjelman alustuksia

(Koodailen lähinnä itsekseni, joten käyttämäni terminologia ei välttämättä ole ihan yleisen käytännön mukaista.)

Käytän jupyterlab nimistä ohjelmointiympäristöä. Ohjelmointiin se on pikkuisen kömpelö, mutta sillä voi kätevästi tehdä tällaisia dokumentteja: Yhdistää tekstiä ja matemaattista esitystä python-ohjelmaan ja sen tulosteisiin. Minun openSUSE linuxiini jupyterlabin sai installoitua vetämällä ohjelmistovalikossa rastin ruutuun.

Ohjelman aluksi esittelen aliohjelmakirjastot, joita haluan käyttää:

import matplotlib.pyplot as plt
import matplotlib.patches as patch
#from matplotlib.collections import PatchCollection
%matplotlib inline
import math as math
import numpy as np
import sympy as sp
from IPython.display import display

sp.init_printing()

• matplotlib.pyplot-kirjastossa on kuvaajien piirtelyyn sopivia funktioita. import-komennon jälkeen voin kutsua esimerkiksi pyplot-kirjaston funktiota plot seuraavasti: plt.plot(x,y)
• matplotlib.patches kirjastossa on kaaviokuvien piirtoon sopivia funktioita.
• %matplotlib inline komento kertoo, että haluan ohjelmani piirtämät kuvaajat näkyviin tähän samaan dokumenttiin.
• math-kirjaston funktiota voi kutsua seuraavasti: y = math.sin(pi)
• numpy-kirjastossa on erilaisia numeerisen laskennan funktiota ja sympy-kirjastossa symbolisen laskennan funktiota.

init_printing() liittyy jotenkin symbolisten lausekkeiden esittämiseen. Olen joskus siihen perehtynyt ja hyödylliseksi todennut.

Integrointirutiineja

Täällä on johdantoa differentiaaliyhtälöiden numeeriseen ratkaisemiseen

Tuo juttu antoi käsillä olevaa tehtävää varten riittävän, 'kansantajuisen' selityksen differentiaaliyhtälöiden ratkaisemisesta. Toivottavasti se herätti kiinnostuksen oikeaan matemaattiseen esitykseen siitä, mitä ovat derivaatta, integraali ja differentiaaliyhtälöt. On tärkeää oppia ymmärtämään oikea matemaattinen tapa osoittaa ratkaisu oikeaksi.

Differentiaaliyhtälöiden ratkaisemiseen on yleensä järkevää käyttää matematiikkakirjastoissa olevia valmiita rutiineja, mutta koodasin tähän uteliaisuudesta ja harjoituksen vuoksi yhden niistä, Runge-Kutta algoritmin

yy(t+dt) = rk4(fdydt,yy,uu,dt)

• fdydt(yy,uu): käyttäjän kirjoittama funktio, joka palauttaa tilamuuttujien yy derivaatat tilan yy ja ulkoisten ohjausten uu funktiona.
• yy: järjestelmän tilamuuttujien $\boldsymbol{y}\,(t)$ lista (=taulukko, vektori, ...)
• uu: $\boldsymbol{u}\,(t)$: ulkoiset ohjaukset ja parametrit
• dt: integrointiaskel
• funktio rk4 palauttaa listan yy1, joka on $ \boldsymbol{y}\,(t+dt)$

Kirjoittamallani rutiinilla voi siis ratkaista usean muuttujan differentiaaliyhtälöryhmiä. Myöhemmin näemme, miten se tapahtuu.

def rk4(fdydt, yy, uu, dt):
    kk1 = fdydt(yy, uu)
    yk1 = [yy[i] + k1 * dt / 2.0 for i, k1 in enumerate(kk1)]
    kk2 = fdydt(yk1, uu)
    yk2 = [yy[i] + k2 * dt / 2.0 for i, k2 in enumerate(kk2)]
    kk3 = fdydt(yk2, uu)
    yk3 = [yy[i] + k3 * dt for i, k3 in enumerate(kk3)]
    kk4 = fdydt(yk3, uu)
    yy1 = [y + dt / 6.0 * (kk1[i] + 2.0 * kk2[i] + 2.0 * kk3[i] + kk4[i]) for i, y in enumerate(yy)]
    return yy1

Metsän kasvu, hiilikierto, lämpeneminen

Kaaviokuva

Piirretään kaaviokuva ilmiöiden välisistä riippuvuuksista ja kirjoitetaan riippuvuudet sen jälkeen differentiaaliyhtälöiksi.

Jos on harjoitellut piirtämistä, havainnollisimman ja mieleenpainuvimman kaaviokuvan saanee vapaalla kädellä piirtämällä. Teknisesti täsmällisen näköisen, mutta 'kuivakkaamman' kuvan saa kaavioiden piirtelyyn tarkoitetuilla ohjelmilla, vaikkei ihan taiteilija olisikaan.

Tässä kokeilin kolmatta vaihtoehtoa: Kaavion piirtämistä ohjelmallisesti. Yksittäisen kaavion piirtämiseen — ainakin ensimmäisellä kerralla — tämä on hitain ja kömpelöin tapa. Hyödyllinen tämä on, jos joutuu piirtelemään paljon samantyylisiä kaavioita varsinkin jos kaaviossa ilmaistava data tuotetaan ohjelmallisesti.

Seuraavaa koodin pätkää ei kannata lukea, ellei aio itse kokeilla jotain vastaavaa. Siinä tapauksessa kannattaa hakukonettaa 'matplotlib.patches'

# Kaavion parametrit olisi ollut hienoa laittaa funktion parametreiksi tai tehdä kaavioista olio

def fkaavio():
    # Vasemmalla alhaalla olevan laatikon vasen alakulma
    x0 = 1.0 
    y0 = 1.0
    # alarivin laatikoiden leveys ja korkeus
    w0 = 7.0
    h = 2.5
    # alaririn laatikoiden laatikoiden välimatka
    dx = 6.0
    # Ylälaatikon vasen alakulma
    x1 = x0 + dx
    y1 = y0 + 4.0*h
    # ylälaatikon leveys
    w1 = 3.5*w0 + 3.0*dx
    
    # nuolten koordinaatit olisi ollut tyylikästä laittaa nekin tähän

    # perusväreillä on nimet, mutta värit voi antaa myös RGB muodossa (red, green, blue)
    # missä red, green ja blue ovat väliltä [0.0, 1.0]
    
    fig, kaavio = plt.subplots(figsize=(14, 8), facecolor=(0.9, 1.0, 1.0)) 

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0, y0+0.5*h), w0, 0.3*h, boxstyle='round', facecolor='r', edgecolor='black'))    
    kaavio.text(x0+0.25*w0, y0+0.6*h, "T1", fontsize='x-large')

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0, y0-0.75*h), w0, h, boxstyle='round', facecolor='r', edgecolor='black'))    
    kaavio.text(x0+0.25*w0, y0-h/2, "T2", fontsize='x-large')

    kaavio.add_patch(patch.Ellipse(
        (x0+w0/2, y0+h/4), 2.25*w0, 3.5*h, 
         edgecolor='r', facecolor='none', linewidth=2))
    
    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+3.0*w0/4.0, y0+h), posB=(x0+3.0*w0/4.0, y1+3.0*h),
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='r'))
    kaavio.text(x0+0.1*w0, y1+3.5*h, "q_in", fontsize='x-large')

    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+w0/4.0, y1+3.0*h), posB=(x0+w0/4.0, y0+h), 
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='r'))
    kaavio.text(x0+0.6*w0, y1+3.5*h, "q_out", fontsize='x-large')

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0+w0/2.0, y1), w1, h, boxstyle='round', color='b', alpha=0.5))
    kaavio.text(x1+0.3*w1, y1+0.4*h, "ilma, CO2", fontsize='x-large')

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0+w0+dx, y0), w0, h, boxstyle='round', color='g'))
    kaavio.text(x0+w0+dx+0.05*w0, y0+0.4*h, "metsä CO2", fontsize='x-large')
    
    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+1.5*w0+dx, y1), posB=(x0+1.5*w0+dx, y0+h),
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='b'))
    kaavio.text(x0+1.25*w0+1.0*dx, y0+2.5*h, 
                "kasvu", fontsize='x-large',bbox=dict(facecolor='w', alpha=0.5, linewidth=0))

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0+2*w0+2.0*dx, y0), w0, h, boxstyle='round', color='y'))
    kaavio.text(x0+2*w0+2.0*dx,  y0+0.4*h, "varasto CO2", fontsize='x-large')
    
    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+2.0*w0+dx, y0+0.5*h), posB=(x0+2.0*w0+2.0*dx, y0+0.5*h), 
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='y'))
    kaavio.text(x0+2.0*w0+1.1*dx, y0+0.6*h, "hakkuut", fontsize='x-large')

    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+2.5*w0+2.0*dx, y0+h), posB=(x0+2.5*w0+2.0*dx, y1), 
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='b'))
    kaavio.text(x0+2.25*w0+2.0*dx, y0+2.5*h, "poltto", 
                fontsize='x-large',bbox=dict(facecolor='w', alpha=0.5, linewidth=0))

    kaavio.add_patch(patch.FancyBboxPatch(
        (x0+3*w0+3.0*dx, y0-2*h), w0, h, boxstyle='round', color='black'))
    kaavio.text(x0+3*w0+3.0*dx+0.1*w0,  y0-1.6*h, "fossiiliset", fontsize='x-large', color='w')

    kaavio.add_patch(patch.FancyArrowPatch(
        posA=(x0+3.5*w0+3.0*dx, y0 - h), posB=(x0+3.5*w0+3.0*dx, y1), 
        arrowstyle="fancy, head_length=40.0, head_width=40.0, tail_width=20.0", color='b'))
    kaavio.text(x0+3.4*w0+3.0*dx, y0+2.0*h, "fos", 
                fontsize='x-large',bbox=dict(facecolor='w', alpha=0.5, linewidth=0))


    #Use adjustable='box-forced' to make the plot area square-shaped as well.
    kaavio.set_aspect('equal', adjustable='datalim')
    plt.axis('off')
    kaavio.plot()   #Causes an autoscale update.
    plt.show()
    plt.close()

fkaavio()

Kaavio esittää lämmön ja hiilen kiertoa maapalloa muistuttavalla planeetalla. Energialähteinä käytetään puuta ja fosiilisia polttoaineita. Simuloinnissa tarkasteltava tehon tarve voidaan sovittaa halutuksi prosenttiosuudeksi puun maksimikasvusta, koska näin saadaan parhaiten esille se, mitä merkitystä on sillä, että osa fossiilista polttoaineista korvataan puulla. Simuloinnin loppupuolella tehon tarve lasketaan nollaan, eli oletetaan tuulen, auringon, säästämisen yms. vähentävän polttamisen tarvetta.

Tarkastelun alussa planeetan metsät ovat vastaistutettua taimikkoa.

• Aurinko lämmittää ($q_{in}$) ilmaa ja planeetan pintaa ($T_1$). Osa lämmöstä siirtyy hitaasti syvemmälle maahan ja meriin ($T_2$);
• Planeetta jäähtyy säteilemällä ($q_{out}$). Jäähdytys riippuu ilman ja planeetan pinnan lämpötilasta sekä hiilidioksidin määrästä ilmakehässä. Hiilidioksidi sitoo itseensä osan planeetan säteilemästä;
• Metsän kasvu sitoo hiilidioksidia ilmakehästä;
• Metsästä korjattu puu siirtyy varastoon. Puun polttaminen varastosta tuottaa CO2:ta ilmakehään;
• Fosiilisten polttoaineiden polttaminen tuottaa CO2:ta ilmakehään;
• Ainakin toistaiseksi maapallolla ilmakehän hiilidioksidista osa liukenee meriveteen.

Tämä tarkastelussa ei oteta huomioon metsien virkistysarvoja eikä pohdita metsien kykyä sopeutua muuttuvaan ilmastoon.

Jos haluttaisiin tarkemmin vertailla avohakkuuta ja jatkuvaa kasvatusta, mallin puut eivät saisi olla yhtä massaa vaan yksilöitä, joiden yksi ominaisuus on ikä. Sellaistakin olen kokeillut

Sama yhtälöinä

$$kasvu(t) = c_{kasvu} \cdot \left( CO2_{kypsa} - CO2_{metsa}(t) \right) \cdot CO2_{metsa}(t)$$

$kasvu(t)$ kuvaa käsitystäni metsään puiden kasvun myötä sitoutuvasta hiilidioksidista. Arvelen nuoren yksin kasvavan puun (ja siihen sitoutuneen hiilen) kasvuvauhdin verrannolliseksi puun kokoon, koska se on verrannollinen yhteyttävän lehtivihreän määrään sekä ravinteita ja vettä sitovan juuriston määrään.

Metsän tihetessä puut alkavat varjostaa toisiaan ja kilpailla vedestä ja ravinteista. Kasvuvauhti hidastuu. Yksittäisen puunkin kasvuvauhti taitaa vanhetessa hidastua, vaikka saisi yksinään kasvaa. Laitoin kasvuvauhdin verrannolliseksi erotukseen 'täysikasvuisen' metsän hiilimäärästä. Näin metsän kasvu hidastuu olemattomaksi ylärajaa lähestyttäessä. En tiedä juuri mitään metsän ja puiden kasvun biologiasta, joten mallini on sangen huteralla pohjalla.

Joidenkin artikkelien mukaan puolet metsien sitoutuneesta hiilestä on sitoutunut maaperään. Lahottajabakteerit tuottavat hiilidioksidin lisäksi hiiliyhdisteitä, jotka jäävät maaperän humuskerrokseen. Monissa artikkeleissa on esitetty, että avohakkuun jäljiltä metsämaa on aluksi hiilen lähde. Tästä en kuitenkaan varsinaisesti tiedä mitään enkä sitä sisällyttänyt malliini.

Metsän sitoutuneen hiilidioksidimäärän muutosnopeus $$\frac{d}{dt}metsa_{CO2}(t) = kasvu(t) - hakkuu(t)$$

Metsästä korjattuun puuhun sitoutuneen hiilidioksimäärän muutosnopeus $$\frac{d}{dt}varasto(t) = hakkuu(t) - poltto(t)$$

Aina kun metsä saavuttaa hakkuukypsyyden — puun määrä kasvaa yli asetetun arvon — siitä hakataan varastoon tietty prosenttiosuus.

Esimerkiksi rakentamiseen käytetty puu on varastoitu pitkäaikaisesti.

Hiilivarastoa kuluttaa polttaminen. Pehmopaperista iso osa päätynee nopeahkosti energiajakeeksi. Polttohakkeen ja biodieselin hiili pöläytetään taivaalle melko nopeasti hakkuun jälkeen. Tässä ohjelmassa varastoa poltetaan vakiovauhdilla, mikäli poltettavaa on, eli kun $varasto(t) > 0$, energiatarve tyydytetään puulla, muuten poltetaan fosiilisia.

Tehon tarve sovitetaan halutuksi prosenttiosuudeksi metsän maksimikasvusta. Simuloinnin loppupuolella energian tarve alkaa laskea nollaan. Ehkä aletaan saada uusiutuvia käyttöön ja on löydetty säästökeinoja.

Ilmakehässä olevan hiilidioksidimäärän muutosnopeus $$\frac{d}{dt}CO2_{ilma}(t) = poltto(t) - kasvu(t) + fos(t)$$

Planeetalta avaruuteen säteilevä lämpö $$q_{out}(t) = c_{jaahdytys}T_1(t) - c_{greenh} CO2_{ilma}(t)$$

Säteily riippuu planeetan pintakerrosten lämpötilasta. Lineaarinen riippuvuus ei vastaa todellisuutta, mutta on sentään oikeansuuntainen. Maapallolla säteily heikkenee ilmaston lämmetessä, koska säteilyä hyvin heijastavan jään ja lumen pinta-ala vähenee. Muitakin vaikutuksia suuntaan ja toiseen ilmaston lämpenemisellä on.

Tästä säteilystä sitä suurempi osa jää ilmakehään, mitä enemmän ilmakehässä on hiilidioksidia.

Pinnasta syvemmälle maaperään ja meriin siirtyvä lämpö $$q_{12}(t) = c_{12} \left( T_1(t)- T_2(t) \right)$$

Lämmön varastoituminen $$cM_1 \frac{d}{dt}T_1(t) = q_{in}(t) - q_{12}(t) - q_{out}(t)$$ $$cM_2 \frac{d}{dt}T_2(t) = q_{12}(t)$$

Koska simuloinnissa tarkastellaan poikkeamaa tasapainotilasta, planeetalle keskustähdestä lankeava säteily on simuloinnissa asetettu $q_{in}(t) = 0$

$c M_1 T_1(t)$ on ilmakehään ja planeetan pintakerroksiin varastoitunut lämpö. $c M_2 T_2(t)$ on syvemmälle maahan ja meriin varastoitunut lämpö.

En käytä suureiden todellisia arvoja vaan skaalaan muuttujien arvot niin, että 0 vastaa sitä, millaista oli "vanhaan hyvään aikaan" ja 100:n voi ajatella olevan siedettävän yläraja.

Olisi ehkä tyylikkäämpää ja selkeämpää, jos lisäisin yhtälöihin kertoimet niin, että puun määrän mittana olisi kiloja tai kuutioita, hiilidioksin määrä kiloja, lämpötilat asteita ja muutenkin parametrit kohdallaan. Hyvä lisäharjoitus.

Mallini ei kuvaa todellisuutta kovinkaan tarkasti, vaan lieneekö kaikilla webiin kirjoittajilla mielessään näinkään monipuolista ja tarkkaa mallia, jonka perusteella he väitteitään esittävät.

Simulointimallin koodi

Tämä ohjelma on kasvanut pienestä prototyypistä kymmenien iteraatioiden kautta nykyiselleen — vähän isommaksi prototyypiksi. Valmis tämä on sitten, kun kyllästyn jatkamaan.

Tässä ohjelmassa lopussa tehdään monta simulaatiota, joiden koodit ovat melkein toistensa kopioita. Olisi mukava tehdä useita simulaatioita pikkuisen eri parametreillä ja piirtää tulokset päällekkäin, että näkisi, kuinka herkkä malli on parametreille. Tätä varten pitäisi kuitenkin tehdä funktio simulaatio(), jolle annettaisiin mallin parametrit argumenttina. Muutakin paranneltavaa on paljon.

Mutta toistaiseksi ohjelma on siis valmis tällaisenaan ;-)

class Simu:
    def __init__(self, Tf):
        self.Tf = Tf      # simuloinnin loppuaika 'vuosina' 
        self.dt = 1.0/12.0  # simuloinnissa käytettävä integrointiaskel. 1/12 = kuukausi.(tarpeettoman lyhyt)
        self.dim = int(self.Tf/self.dt)  # aikapisteiden lukumäärä

Metsän kasvu ja hakkuut

Tarkistetaan määrävuosina, onko metsä hakkuukypsää ja jos on, hakataan puista tietty osuus.

Käytännössä siihen, minkä ikäiset puut kannattaa hakata vaikuttanee metsän tuottovaade: "Kannattaako puut hakata nyt vai 10 vuoden kuluttua, jos tuottovaade on 4%? ja puiden arvon arvioidaan kasvavan 10 vuodessa 1.4 kertaiseksi?"

# Korkoa korolle, kymmenen vuotta, neljän prosentin tuottovaade
print('Metsän arvon pitäisi kasvaa ',(1+0.04)**10, ' kertaiseksi')

Metsän arvon pitäisi kasvaa  1.4802442849183444  kertaiseksi

Ei kannattane jättää kasvamaan. Talouslaskennasta tiedän niin vähän, etten yritä tämän enempää.

Allaolevaa funktiota hakkuukutsutaan ajan dt-välein eli joka simulointiaskeleella. Siksi se saattaa näyttää omituiselta.

Mikäli metsän puumäärä m on yli hakkuurajan, aloitetaan hakkuut. Hakkuut lopetetaan, kun hakkuukesto simulointiaskelta on kulunut. Jos simulointiaskel on kuukausi, hakkuu siis jatkuu hakkuukesto kuukautta.

class Metsa:
    def __init__(self):
        self.c_kasvu = 0.01*simu.dt
#        print("nuoren väljästi kasvavan puun kasvuvauhti: ", self.c_kasvu/self.dt*100.0, "prosenttia vuodessa")
        self.co2_kypsa = 100.0   # Metsän maksimikoko
        self.hakkuuprosentti = 99.0  # kun hakataan, hakataan suunnilleen näin monta prosenttia metsästä
        self.hakkuuraja = 70.0   # Käynnistetään hakkuut, jos metsän koko on tämän rajan yli
        self.hakkuukk = 0
        self.maara_kk = 0
        self.hakkuuKesto = 6
       
    def hakkuu(self, m):
        if self.hakkuukk == 0 and m > self.hakkuuraja:
            self.hakkuukk = self.hakkuuKesto
            self.maara_kk = 0.01*self.hakkuuprosentti*m*12/self.hakkuuKesto
            # print('hakkuu alkaa! m: ', m)
        if self.hakkuukk > 0:
            # print('hakkuukk: ', self.hakkuukk, 'm: ', m, 'hakkuumäärä: ', self.maara_kk)
            self.hakkuukk -= 1
            return self.maara_kk
        return 0.0

    def kasvu(self, co2_metsa):
        kasvu = self.c_kasvu*(self.co2_kypsa - co2_metsa)*co2_metsa
        return kasvu

Poltto

Yritetään korvata fossiilisia polttoaineita puulla. Simuloinnin loppupuolella energian tarpeen oletetaan vähenevän säästämisen sekä aurinko- ja tuulivoiman lisääntymisen myötä niin, että hakattua puuta voidaan käyttää vaikka rakentamiseen.

Energian tuotannon yhtälöissä ei ole otettu huomioon kummankaan polttoaineen tuotannossa syntyviä CO2-päästöjä. Näin summittaisissa tarkasteluissa sellaisella ei ehkä olisi suurtakaan merkitystä.

class PuuVarasto:
    def __init__(self):
        self.kulutus = 0.75*metsa.kasvu(50.0)
        print("Energian kulutus alussa: ", self.kulutus)
        self.bioOsuus = 1.0  # Puun osuus polttoaineesta, kun puuta on varastossa
        self.t1 = 0.6*simu.Tf # Tästä hetkestä eteenpäin energian tarve alkaa laskea.
        self.t2 = self.t1 + 30.0 # Energian tarve laskee nollaan hetkeen t2 mennessä

    def poltto(self, vuosi, puuVarasto):
        if vuosi < self.t1:
            kulutus = self.kulutus
        elif vuosi > self.t2:
            kulutus = 0.0
        else:
            kulutus = self.kulutus/(self.t1 - self.t2)*vuosi - self.kulutus*self.t2/((self.t1-self.t2))
            
        if puuVarasto > 0:
            bio = kulutus*self.bioOsuus
            fos = kulutus - bio
            return (bio, fos)
        return (0.0, kulutus)

Planeetan lämpeneminen

Yhtälöt uudelleen: $$q_{out}(t) = c_{jaahdytys}T_1(t) - c_{greenh} CO2_{ilma}(t)$$ $$q_{12}(t) = c_{12} \left( T_1(t)- T_2(t) \right)$$ $$cM_1 \frac{d}{dt}T_1(t) = q_{in}(t) - q_{12}(t) - q_{out}(t)$$ $$cM_2 \frac{d}{dt}T_2(t) = q_{12}(t)$$

class Maapallo:
    def __init__(self):
        self.m1_pallo = 20.0      # Lämpökapasiteetti
        self.m2_pallo = 100.0      # Lämpökapasiteetti
        self.c_jaahdytys = 0.8
        self.c_greenh = 1.2
        self.c_12 = 1.5
        self.q_in = 0.0

    def lammitys(self, T1, T2, co2_ilma):
        q12 = self.c_12*(T1 - T2)
        q_out = self.c_jaahdytys*T1 - self.c_greenh*co2_ilma
        dT1 = (self.q_in - q12 - q_out)/self.m1_pallo 
        dT2 = q12/self.m2_pallo 
        return (dT1, dT2)

Systeemin tilayhtälö

Integrointirutiinien vaatima järjestelmän tilamuuttujien aikaderivaatat palauttava funktio. Yhtälöt uudelleen:

$$ kasvu(t) = c_{kasvu} \cdot \left( CO2_{kypsa} - CO2_{metsa}(t) \right) \cdot CO2_{metsa}(t) \\ \frac{d}{dt}metsa_{CO2}(t) = kasvu(t) - hakkuu(t) \\ \frac{d}{dt}varasto(t) = hakkuu(t) - poltto(t) \\ \frac{d}{dt}CO2_{ilma}(t) = poltto(t) - kasvu(t) + fos(t) \\ q_{out}(t) = c_{jaahdytys}T_1(t) - c_{greenh} CO2_{ilma}(t) \\ q_{12}(t) = c_{12} \left( T_1(t)- T_2(t) \right) \\ cM_1 \frac{d}{dt}T_1(t) = q_{in}(t) - q_{12}(t) - q_{out}(t) \\ cM_2 \frac{d}{dt}T_2(t) = q_{12}(t) \\ $$

Tilamuuttujia, joiden derivaatat lasketaan, ovat siis: $CO2_{metsa}(t)$, $CO2_{ilma}(t)$, $T_1(t)$, $T_2(t)$ ja $varasto(t)$. Järjestelmään vaikuttavat ulkopuoliset 'ohjaukset' ovat hakkuut ja puun ja fossiilisten poltto.

def fdxdt(xx, uu):
    [co2_metsa, co2_ilma, T1, T2, puuVarasto] = xx
    [hakkuut, bio, fos] = uu
    kasvu = metsa.kasvu(co2_metsa)
    dco2_metsa = kasvu - hakkuut
    dco2_ilma = bio - kasvu + fos
    (dT1, dT2) = maapallo.lammitys(T1, T2, co2_ilma)
    dVarasto = hakkuut - bio
    return [dco2_metsa, dco2_ilma, dT1, dT2, dVarasto]

Tulosten piirtäminen

Ketä kiinnostaa, voi alta selvittää, miten piirtää kuvaajia pythonilla.

# CO2, uu ja yy ovat 2-dimensioisia taulukoita esitettynä listojen listoina
def plottaa(tt, CO2, uu, yy, otsikko):
    nyky_alku = 30
    nyky_loppu = 120
    # Matplotlib-kirjaston vaatimuksesta listoista pitää tehdä numpy-taulukoita
    CO2 = np.array(CO2)
    uu = np.array(uu)
    yy = np.array(yy)
    tt = np.array(tt)

    # Toiseen koordinaatistoon tilamuuttujat, toiseen ohjaukset ja kasvu
    # Neljä koordinaatistoa
    # 0,0 tilamuuttujat koko simuloinnin ajalta
    # 0,1 ohjaukset ja 'mitatut' suureet koko simuloinnin ajalta
    # 1,0 tilamuuttujat 'lähimenneisyydestä ja -tulevaisuudesta'
    # 1,1 ohjaukset ja 'mitatut' suureet 'lähimenneisyydestä ja -tulevaisuudesta'
    f, axarr = plt.subplots(2, 2, figsize=(16, 12))

    plt.suptitle(otsikko, color='black')
    axarr[0,0].set_ylim([-25.0, 125.0])
    axarr[0,0].plot(tt, CO2[:, 0], 'g-', label='metsä CO2')
    axarr[0,0].plot(tt, CO2[:, 1], 'b-', label='ilma CO2')
    axarr[0,0].plot(tt, CO2[:, 2], 'r', linestyle='--', label='T1')
    axarr[0,0].plot(tt, CO2[:, 3], 'r', linestyle=':', label='T2')
    axarr[0,0].plot(tt, CO2[:, 4], 'k', linestyle=':', label='puuvarasto CO2')
    axarr[0,0].grid(True)
    axarr[0,0].tick_params(axis='both', colors='black')
    axarr[0,0].legend(loc='upper left')
    
    axarr[0,1].plot(tt, uu[:, 1], 'r', linestyle='--', label='poltto')
    axarr[0,1].plot(tt, uu[:, 2], 'black', linestyle=':', label='fosiiliset')
    axarr[0,1].plot(tt, yy[:, 0], 'g', label='kasvu')
    axarr[0,1].grid(True)
    axarr[0,1].tick_params(axis='both', colors='black')
    axarr[0,1].legend(loc='upper right')

#    print("\n\n  Lähimenneisyys ja -tulevaisuus")

    ymin = -1.0
    ymax = 1.0
    for i in range((nyky_alku-1)*12, (nyky_loppu+1)*12):
        for j in range(0,5):
            ymin = min(ymin, CO2[i][j])
            ymax = max(ymax, CO2[i][j])
        
    axarr[1,0].set_ylim([ymin, ymax])
    axarr[1,0].set_xlim([nyky_alku, nyky_loppu])
    axarr[1,0].plot(tt, CO2[:, 0], 'g-', label='metsä CO2')
    axarr[1,0].plot(tt, CO2[:, 1], 'b-', label='ilma CO2')
    axarr[1,0].plot(tt, CO2[:, 2], 'r', linestyle='--', label='T1')
    axarr[1,0].plot(tt, CO2[:, 3], 'r', linestyle=':', label='T2')
    axarr[1,0].plot(tt, CO2[:, 4], 'k', linestyle=':', label='puuvarasto CO2')
    axarr[1,0].grid(True)
    axarr[1,0].tick_params(axis='both', colors='black')
    axarr[1,0].legend(loc='upper left')
    
    ymin = -1.0
    ymax = 1.0
    for i in range((nyky_alku-1)*12, (nyky_loppu+1)*12):
        for j in range(1,3):
            ymin = min(ymin, uu[i][j])
            ymax = max(ymax, uu[i][j])
        ymin = min(ymin, yy[i][0])
        ymax = max(ymax, yy[i][0])

    axarr[1,1].set_ylim([ymin, ymax])
    axarr[1,1].set_xlim([nyky_alku, nyky_loppu])
    axarr[1,1].plot(tt, uu[:, 1], 'r', linestyle='--', label='poltto')
    axarr[1,1].plot(tt, uu[:, 2], 'black', linestyle=':', label='fosiiliset')
    axarr[1,1].plot(tt, yy[:, 0], 'g', label='kasvu')
    axarr[1,1].grid(True)
    axarr[1,1].tick_params(axis='both', colors='black')
    axarr[1,1].legend(loc='upper right')
    
    plt.show()
    plt.close()

Simuloinnit

Alkuarvoja

# Luodaan oliot
simu = Simu(200)
metsa = Metsa()
maapallo = Maapallo()
varasto = PuuVarasto()

Energian kulutus alussa:  1.5624999999999998

metsa_0 = 0.2   # Metsän CO2-ekvivalentti puumäärä CO2 simuloinnin alussa
ilma_0 = 0.0  # Ilman CO2-määrä simuloinnin alussa
T1_0 = 0.0    # lämpötila simuloinnin alussa
T2_0 = 0.0    # lämpötila simuloinnin alussa
bio_0 = 0.0  # puun varasto alussa

# Luodaan luettelo — python-kielellä lista — aikapisteistä, joissa lasketaan muuttujien arvot
print("aikapisteiden lukumäärä: ", simu.dim)
tt = [i * simu.dt for i in range(0, simu.dim)]  # aikapisteet

# Luodaan listat, joihin talletetaan tilamuuttujien ja ohjausten arvot kussakin aikapisteessä.
CO2 = [[0,0,0,0,0] for i in range(0, simu.dim)]  # Tilamuuttujat
uu = [[0, 0, 0] for i in range(0, simu.dim)]  # Ohjaukset
yy = [[0] for i in range(0, simu.dim)]  # Muut piirrettävät suureet

aikapisteiden lukumäärä:  2400

Simulaatio: CO2 vaikutus lämpötilaan

Mitä tapahtuu, kun CO2 ensin kasvaa, sitten käännetään laskuun? Puut eivät kasva, eikä niitä polteta.

CO2[0] = [metsa_0, 0.0, 0.0, 0.0, 0.0]  # tilamuuttujien arvot alkuhetkellä t = 0

# Nollataan kasvu tämän simuloinnin ajaksi
Tmp1 = metsa.c_kasvu
metsa.c_kasvu = 0.0
kaanne1 = 600
kaanne2 = 1400
kaanne3 = 2000
# Lasketaan tila(i+1) tila(i):n ja ohjausten(i) avulla
# CO2 pakotetaan tiettyyn arvoon kullakin hetkellä
for i in range(0, simu.dim - 1):
    if i < kaanne1:
        CO2[i][1] = i*100/kaanne1
    elif i > kaanne2:
            CO2[i][1] = max(0, 100.0 - 100.0*(i-kaanne2)/(kaanne3-kaanne2))
    else:
        CO2[i][1] = 100.0
    # CO2[i][1] = 100.0  # testi, poista
    uu[i] = [0.0, 0.0, 0.0]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])
plottaa(tt, CO2, uu, yy, "puu ei kasva")

metsa.c_kasvu = Tmp1

Tulos vastaa jotenkuten odotuksiani, ehkä. Vaikka kasvihuonekaasujen nettopäästöt saataisiin nollaan, maapallon lämpenisi vielä pitkään, koska se ei ole ehtinyt lämmentä samaa vauhtia kuin kasvihuonekaasut ilmakehässä ovat lisääntyneet. Jäähtyminen alkaa tietysti vasta, jos kasvihuonekaasujen pitoisuus ilmakehässä onnistutaan vähentämään.

Todellisuudessa maapallon lämpeneminen ja lämmön varastoituminen ovat paljon malliani monimutkaisempia ilmiöitä. Pahimmassa tapauksessa napajäätiköiden sulaminen ja muut säteilyn heijastumista vähentävät ilmiöt saattavat aiheuttaa sen, että 'ylimääräisten' kasvihuonekaasujen imaiseminen pois ei riitä palauttamaan lämpötilaa ennalleen.

Tehtävä: Tarkista IPCC:n raportista, millä dynamiikalla lämpeneminen oikeasti tapahtuu.

Simulaatio: Annetaan metsän kasvaa vapaasti

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, bio_0]  # tilamuuttujien arvot alkuhetkellä t = 0
hakkuu = 0.0
varasto.bioOsuus = 0.0

for i in range(0, simu.dim - 1):
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi, CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])

plottaa(tt, CO2, uu, yy, "puuta ei hakata eikä polteta")

varasto.bioOsuus = 1.0

Metsää ei hakata ja energia tuotetaan fosiilisilla. Planeetan lämpötila karkaa tietysti käsistä.

Nykytilanteessa, välillä 30-80, metsät ovat hyvässä kasvussa ja auttavat pitämään ilmastonmuutosta kurissa.

Hakattu puu poltetaan

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, 0.0]  # tilamuuttujien arvot alkuhetkellä t = 0
for i in range(0, simu.dim - 1):
    hakkuu = metsa.hakkuu(CO2[i][0])
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi, CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])
plottaa(tt, CO2,uu, yy, "puu poltetaan, avohakkuu")

Avohakkuu ja puun poltto on lähitulevaisuudessa haitaksi, koska avohakkuu pysäyttää metsän kasvun, mutta siitä eteenpäin puun polttaminen on tietysti fossiilisten polttamista parempi vaihtoehto. Lämmin tulee kuitenkin.

Koska mallissa kasvuvauhti riippuu puun määrästä, malli tuottaa hakkuiden alussa epärealistisen kasvupiikin.

Kaikkia planeetan metsiä ei tietenkään avohakattaisi samanaikaisesti vaan palsta kerrallaan, jolloin avohakkuun aiheuttama kasvun pysähdys jakautuisi eri vuosille. Haitta olisi yhtä suuri, mutta peittyisi keskiarvon alle.

hakattu puu varastoidaan

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, bio_0]  # tilamuuttujien arvot alkuhetkellä t = 0

varasto.bioOsuus = 0.0

for i in range(0, simu.dim - 1):
    hakkuu = metsa.hakkuu(CO2[i][0])
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi,CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])
plottaa(tt, CO2, uu, yy, "kaikki hakattu puu varastoidaan, avohakkuu")

varasto.bioOsuus = 1.0

Sama tulos kuin yllä. Tätä jouduin hetken hämmästelemään. Näyttää yhdentekevältä, poltetaanko hakattu puu vai varastoidaanko puu ja poltetaan sen sijasta fossiilisia, kun ei oteta huomioon polttoaineiden tuotannon CO2-päästöjä. Loogista loppujen lopuksi.

Osa poltetaan, osa varastoidaan

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, 0]  # tilamuuttujien arvot alkuhetkellä t = 0

varasto.bioOsuus = 0.5

for i in range(0, simu.dim - 1):
    hakkuu = metsa.hakkuu(CO2[i][0])
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi,CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])
plottaa(tt, CO2,uu, yy, "puuta poltetaan ja varastoidaan")

Sama tulos tästäkin tietysti.

Varovainen hakkuu

Osa puusta poltetaan, osa varastoidaan

metsa.hakkuuraja = 85.0
metsa.hakkuuprosentti = 50.0

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, bio_0]  # tilamuuttujien arvot alkuhetkellä t = 0

for i in range(0, simu.dim - 1):
    hakkuu = metsa.hakkuu(CO2[i][0])
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi, CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])
plottaa(tt, CO2, uu, yy, "osa puusta poltetaan, hakataan kypsää metsää varovaisesti")

Laitetaan hakkuuraja ylemmäksi ja annetaan metsän järeytyä. Tehdään vain poimintahakkuita, hakataan kerralla vain puolet metsästä — vanhimmat ja huonokuntoiset puut. Lämpeneminen pysyy paremmin aisoissa, koska metsässä on keskimäärin enemmän hiiltä varastossa ja kasvu keskimäärin suurempaa. Hakkuun myötä kasvu kiihtyy, ei pysähdy niin kuin avohakkuussa.

Hakkuita tehdään useammin kuin avohakatessa. Lienee kalliimpaa.

Väitetään, että avohakkaamalla uudistettu metsä kasvaa nopeammin kuin jatkuvalla kasvatuksella. En tiedä, mihin tämä perustuu. Kokeita varten pitää ensin kasvattaa metsään, jossa on jatkuvalle kasvatukselle tyypillisesti monen ikäistä puustoa. Avohakkuulla kasvatetusta homogeenisesta puustosta ei sellaista saane kovinkaan nopeasti.

Varovainen hakkuu, puu poltetaan

metsa.hakkuuraja = 85.0
metsa.hakkuuprosentti = 50.0

varasto.bioOsuus = 1.0

CO2[0] = [metsa_0, ilma_0, T1_0, T2_0, bio_0]  # tilamuuttujien arvot alkuhetkellä t = 0

for i in range(0, simu.dim - 1):
    hakkuu = metsa.hakkuu(CO2[i][0])
    vuosi = i*simu.dt
    (bio, fos) = varasto.poltto(vuosi, CO2[i][4])
    uu[i] = [hakkuu, bio, fos]
    CO2[i + 1] = rk4(fdxdt, CO2[i], uu[i], simu.dt)
    yy[i][0] = metsa.kasvu(CO2[i][0])    
plottaa(tt, CO2, uu, yy, "puu poltetaan, hakataan kypsää metsää varovaisesti")

Sama tulos kuin edellisestä.

Lisäharjoituksia

Jos intoa riittää, ilmiön tarkastelua voi jatkaa loputtomiin. Jos haluaa harjoitella haasteellisempia tehtäviä, voi jatkaa aikaisemman pohjalta.

Optimaalisen hakkuun määrittäminen

Hakkuille voisi määrittää strategian, joka minimoisi planeetan lämpötilan maksimin. Optimointimenetelmästä riippuen yksittäisen suureen maksimin minimointi voi olla vaikeaa. Usein tosielämän kannalta olennaisin kriteeri on optimointitehtävissä hankalin minimoitava/maksimoitava. Hyvän tuloksen ehkä antaisi lämpötilan neliön integraali. Tärkeää olisi keksiä, miten kuvata optimointikriteerissä hakkuiden kiinteät kustannukset. Jokainen hakkuukierros maksanee jotain siitä riippumatta, paljonko puuta hakataan.

Esimerkki optimiohjauksen laskemisesta

Metsän simulointi joukkona yksittäisiä puita

Eri-ikäisten puiden erilaisen kasvun mallintamiseen sopii simulointi, jossa metsän jokainen puu on oma olionsa. Tällaistakin kokeilin ja tulokset olivat mielenkiintoisia, mutta eräiltä osin kovin outoja. Yksittäisen puun kasvun mallini ei taida vastata todellisuutta kovinkaan hyvin. Pitäisi perehtyä tarkemmin puiden kasvuun ja mallintaa metsä huolellisemmin.

Tämä työkirja python-koodina

Muuta ohjelmoinnista