01対象と、「厳密」の定義
ブラウザ上で動作する単一HTMLのアクアリウム・シミュレータである。30cmキューブの水草水槽から600cmのリーフタンクまで5槽が同時に走り、水温・pH・KH・アンモニア・亜硝酸・硝酸・溶存酸素・CO2・DOC・塩分が実時間で更新される。三次元描画は計算コアの出力を読むだけの層であり、体の曲がりや歩脚の接地も同一のモデルから導いている。
本実装で「厳密」と呼ぶ条件を3つ定義した。
1. 収支が閉じること
窒素は生成も消滅もしない。餌として投入された窒素は、魚の体・生物膜・デトリタス・水草・藻類・エビ・換水による排出・脱窒による離脱のいずれかに必ず存在する。この収支を毎ステップ台帳に記録し、テストで相対残差 1e-9 未満を要求する。水量についても同様に 1e-6 L 未満を要求する。この制約は検証手段として有効であり、実装上の誤りの半数はここで検出された(§11)。
2. 文献値と推定値を区別すること
全パラメータに確度 A/B/C を付与している。文献は出典のある実測値、推定は挙動を合わせるために設定した値である。UI上にも根拠を表示する。区別を保存しないと、後から数値の由来を追跡できなくなる。
3. 真値と計器を分離すること
内部状態の真値と、利用者が読み取れる値は別に扱う。試薬(テトラ6in1・API)は段階に丸め、発色待ち時間を持ち、pHペンはドリフトし校正を要する。亜硝酸 0.5 mg/L が試薬では 0 とも 1 とも読める、という状況を再現している。これは演出ではなく、実際の飼育において誤判断が生じる主要因だからである。
02窒素:正値性と保存の同時達成
窒素循環は8プール間の輸送として記述できる。
溶存アンモニア / 亜硝酸 / 硝酸 / 易分解デトリタス / 難分解デトリタス / 生物膜 / 濾過槽トラップ / ゼオライト吸着
dci/dt = Σj Rij − Σj Rji + Piext − Diext
問題となるのは時間刻みである。内部 dt は60秒だが、倍速再生時にはより粗い刻みで進めたい。陽的オイラー法では、硝化速度が大きいときに TAN が負値をとる。負の濃度は物理的に無意味であるうえ、Monod 式の分母に入って発散する。
一方、負値を 0 にクリップすると窒素が消失し、収支が閉じない。
採用したのは修正Patankar–Euler法(Burchard, Deleersnijder & Meister 2003)である。破壊項に加えて生成項も未知数側の比 cjn+1/cjn で重み付けし、線形系を解く。
各輸送項について「出る側の減少」と「入る側の増加」が同一量になるため、厳密に保存し、かつ無条件に正である。8×8のガウス消去を毎ステップ実行するコストは、この保証に対して十分小さい。
下図は同一の硝化連鎖を3手法で解いたものである。Δt を大きくすると差が現れる。
03硝化菌の4集団
一般の解説では硝化菌は一括して扱われることが多いが、少なくとも4集団に分けないと観測される現象が再現できない。
| 集団 | 役割 | Ks (mg N/L) | 性質 |
|---|---|---|---|
| AOB(Nitrosomonas) | NH3→NO2 | 大 | 高アンモニア濃度で有利。立ち上げ期に優占 |
| AOA(アンモニア酸化古細菌) | NH3→NO2 | 0.0019 | 極低濃度に強い。光で阻害される。成熟槽で優占 |
| コマモックス | NH3→NO3(完全酸化) | 小 | 淡水で優占しやすい。海水では検出されない |
| NOB(Nitrospira) | NO2→NO3 | 中 | AOBより Ks が大きい |
分離した結果、いくつかの経験則が明示的な記述なしに再現された。
海水版では文献の状況が異なっていた。Sauder et al. 2024 は海水観賞魚水槽8槽すべてでコマモックスを検出しておらず、AOA優占5槽・AOB優占2槽・同程度1槽と報告している。AOA は 60 µE 程度の光で顕著に阻害され、120 µE 連続照射では停止して暗期にも回復しない(Merbt et al. 2012)。そこで面ごとに「光の当たる割合」を保持している(ライブロック0.35・砂0.6・ガラス0.5・濾過槽0)。リーフ照明下では、ライブロック表面の AOA 活性係数が 0.69 まで低下する。
また AOA はカタラーゼを持たないため、過酸化水素 0.2 µM で停止する(Kim et al. 2016)。藻類対策として過酸化水素を投与するとアンモニアが一時的に上昇する挙動をテスト51で確認している(3 mg/L 投入で4時間後に TAN が1.3倍、30時間で復帰)。
04境界層とアルカリ度
境界層による物質移動の上限
生物膜の硝化速度を Monod 式のみで記述すると、成熟槽の1パス除去率が99%となり、TAN が 0.001 mg/L まで低下する。実測はこれほど低くない。欠落していたのは生物膜表面の拡散境界層である。膜の活性がどれだけ高くても、基質が膜面へ到達する速度に上限がある。
導入後は1パス除去率が40%程度、TAN 0.02–0.04 mg/L で定常となり、実測と整合する範囲に収まった。
アルカリ度の正味消費は 7.14 ではなく 3.57
硝化によるアルカリ度消費として広く引用される「7.14 mg CaCO3 / mg N」は、NH4Cl を外部から投入するフィッシュレスサイクリングの場合にのみ成立する。
魚の排泄・有機物のアンモニア化: R−NH2 + H2O → NH4+ + OH− → +1 eq/mol N
正味 −1 eq/mol N = 3.57 mg CaCO3 / mg N
魚が在籍する水槽ではアンモニアがタンパク質の分解に由来するため、その段階でアルカリ度が 1 eq 回復する。それでも大型オスカーの90cm水槽では1日あたり 0.4 dKH 減少する。東京の水道水(KH 2.5)では1週間持たない。
05水草:CO2は飽和値を下げるが初期勾配を変えない
本実装で最も校正に時間を要した箇所である。
当初、水草の光合成を次のように記述した。光の効果 fL と CO2 の効果 fC を独立に計算して乗じる形である。
fL = 1 − exp(−PAR / Ek) fC = CO2 / (K½ + CO2)
P = pmax · fL · fC · f(T) · f(養分)
この実装では、CO2 添加なしの30cm水槽でアヌビアスとウィローモスが枯死した。実物のアヌビアスは CO2 無添加かつ低光量の条件下でも長期間生存する。
P–I 曲線の文献を再確認したところ、誤りが判明した。CO2 制限は飽和光下の上限(プラトー)を下げるが、初期勾配 α = pmax/Ek は変えない。低光量域では光子束そのものが律速であり、CO2 濃度が低くても光子1個あたりの固定量は変化しないためである。乗算形の実装は、律速でない要因に対して二重に制限を課していた。
P = pmax · fC · ( 1 − exp( −PAR / (Ek · fC) ) ) · f(T) · f(養分) · f(健康)
PAR → 0 のとき P → pmax·PAR/Ek (fC が消去され、初期勾配は不変)
下図で CO2 スライダを下げると、修正後の式では低光量側の傾きが保たれるのに対し、乗算形は原点から曲線全体が潰れる。CO2 無添加の水槽(CO2 はおよそ 1–3 mg/L)でアヌビアスを選ぶと、PAR 30 における純光合成が乗算形では修正後の 1/3 程度まで低下する。実装ではこれに健康度のフィードバック(成長が止まる → 健康度が下がる → さらに成長が落ちる)が重なるため、数十日で枯死に至っていた。
重炭酸(HCO3−)を利用できる種は限られる。バリスネリア・マツモ・ハイグロフィラは利用可能(Sand-Jensen & Gordon 1984 ほか)、アヌビアスやブセファランドラは不可である。「CO2 無添加でも育つ水草」の一部はこの差に起因する。
(a) 健康度の判定に瞬時の PAR を用いていたため、夜間は全個体が光不足と判定され、平均健康度が 0.46 まで低下していた。日平均 PAR を保持し、日照時間で積分した量と比較するよう修正。
(b) 成長速度の上限 rgrmax を1日あたりで定義したまま毎ステップ適用していたため、実際に成長が生じる明期のみで見ると上限が過度に厳しかった。
cap = B·rgrmax/(86400·光周期)·Δt に変更。3件はいずれも「水草が痩せる」という同一の症状を呈するため、切り分けに時間を要した。
06藻類:律速は栄養ではなく付着面積
藻類5種(珪藻・緑スポット・糸状藻・黒髭・ラン藻)を導入したところ、栄養が存在する限り増殖が止まらず、酸素を消費してオスカーの成長を阻害した(成長試験が 25.9 cm で不合格)。
欠落していたのは付着可能面積である。藻類は浮遊しているのではなく、ガラス・葉・流木・底床に付着している。面積あたりの上限(3.5–9 g乾重/m²推定)を種ごとの付着先に応じて与えたところ、試験は合格した(29.9 cm)。
そのうえで、健全な水草が存在すると藻類は抑制される。これはアレロパシー項を記述した結果ではなく、光と栄養が先に消費されるためである。水草なし・強光・過給餌の60cm水槽では60日で乾重 0.97 g(ラン藻優占)、同条件でも密植した場合は1/10以下となる。
07魚:重心まわりの剛体回転ではない
魚が向きを変える際、重心まわりに剛体回転するのではなく、頭部を先に切り込ませ、体幹がその軌跡を追従する。これをルールとして直接記述しないことを目標とした。記述すべきは骨格と流体力学であり、頭部先行は結果として生じるべきである。
骨格
- 剛体の頭蓋(吻端〜後頭顆、全長の 0.25–0.30)推定
- 椎骨チェーン。スズメダイ科26・ニザダイ科22(9+13)は科の標準形質文献、オスカー27・ピラニア36は科の範囲からの推定推定。ピラニアの前4椎はウェーバー器官で癒合し可撓性を持たない
- 尾鰭の受動節2つ(筋による駆動を持たず、流体力で撓む)
- 椎間関節の側方可撓性 flex = 0.25 + 0.75·s1.3(頭→尾で増加)推定
重心
超楕円断面(Γ関数から面積係数を得る)の面積分布を Simpson 則で積分し、局所密度に頭蓋 +12%、鰾(u = 0.28–0.65 で局所断面の最大22%、淡水)を与える。得られた重心位置は吻端から全長比 0.34–0.38 であり、Webb 1978(マス 0.36 TL)および Domenici & Blake 1997(0.35–0.40 TL)と整合する文献。体を曲げた状態でも重心が各節の質量和と一致することをテスト52で検証している。
力
水平面の剛体フレーム動力学として、未知数を「重心の横速度」と「頭蓋フレームのヨー角速度」の2つに限定し、各節について以下を加算する。
- 抵抗力理論(Taylor 1952)の法線抗力 ½ρCnS|vn|vn
- 低速域の線形粘性減衰
- 細長体理論(Lighthill 1971)の付加質量 πρ(d/2)²、およびその運動量が体軸に沿って後方へ輸送される対流項
- 尾鰭の揚力 ½ρCLαUS
- 形状変化にともなう角運動量(体を曲げると反動で尾が逆方向へ振れる)
2×2の半陰的連立方程式を毎サブステップ解く(無条件安定)。
この構成では、屈曲した体は抵抗の異方性により自身の中心線に沿って滑るため、旋回が follow-the-leader となる。下図は同じ 90° 転回を両方式で描いたものである。
08進行波の包絡は全長基準である
魚の遊泳運動学で頻繁に引用される Videler (1993) の横振幅包絡は次の形をとる。
当初これを体長(SL、吻端〜尾柄末端)基準として使用していた。x=1 で A=0.10 となるため「尾柄で 0.10 SL」と解釈していた。
検証で数値が合わず、文献の実測(尾鰭端 0.10 TL)と突き合わせた結果、この式は全長(TL、尾鰭端まで)基準であることが判明した。SL 基準で適用すると尾柄の振幅が過大となり(0.10 に対し正しくは 0.062)、波長も過大になる(1 SL に対し 1.22–1.27 SL)。
修正後の実測値(直進 1.2 BL/s)は、鼻先の横振幅 0.014–0.026 L、頭部のヨー ±2.6–3.1°、重心の横振れ ±0.011–0.014 L、尾鰭端 0.066–0.079 L である。Videler の 0.02 / ±4.6° / 0.10 TL に対し、包絡の形状は一致するが尾端はやや小さい。
09ヌマエビ:25°という制約
ミナミヌマエビとヤマトヌマエビの実装で最も有効だった文献は Arnott, Neil & Ansell 1998(Crangon crangon、体長 11–69 mm)である。テールフリップの高速度撮影から、全屈曲時の頭胸部軸と尾部軸のなす角が 25.0±4.9°、最大伸展 128.2±20.3° という実測が得られている。
これを制約として使用した。腹節間の関節は左右一対の顆をもつ蝶番(dicondylic)であり矢状面1自由度、側屈および捻れはほぼ生じない。8体節(頭胸部+A1–A6+尾節)の関節可動域の総和を 155° に配分すれば、180 − 155 = 25° が導かれる。
運動側も文献のスケーリングに合わせた。テールフリップの屈曲時間は体長回帰(0.34 ms/mm)からミナミヌマエビ(TL 25 mm)で 20–30 ms、再伸展 25–40 ms、2–3連打。逃避方向は刺激から ±63° 以内へは向かわず、最終角 161±55°。歩行は同側が前→後のメタクロナル波(位相 0.33)、対側がほぼ逆相であり、10本の歩脚のうち常に3本以上が接地する(テスト59)。
生活史も実装した。脱皮周期は23℃で成体約3週・稚エビ3–4日、1回の脱皮で甲長が1.13倍になる。抱卵数 21–51個、孵化までの日数は 24/28/32℃ でそれぞれ 21/16/12 日(Tropea 2015)、性比は温度依存(20/23/26℃ で ♀ 75/60/40%)。そしてヤマトヌマエビのゾエアは 8.5–17 ppt 以上の汽水を要求するため、淡水では繁殖しない。この非対称性をモデルに含めると、「ミナミは増えヤマトは増えない」が自然に生じる。60cm水草水槽にミナミ15匹を導入した場合、180日で29–30匹(脱皮106回)となる。
10水流:MAC格子による非圧縮流れ
フィルターによる水流は非圧縮流体として解いている。約7000セルの MAC 格子、半ラグランジュ MacCormack 移流 → 運動量源 → Gauss–Seidel 圧力射影という構成で、1ステップあたり 4 ms(60cm水槽)である。
運動量源は3種類ある。吐出口の噴流(Q·ujet を噴流域の体積で除した加速度)、気泡プルームの浮力(空気流量 × g / (上昇速度 25 cm/s × 断面積))、ヒーターの熱対流(gβΔT ≈ 0.12 cm/s²)である。取水口と吐出口は「速度を規定した固体セル」として扱い、面の流量総和が Q に一致するよう構成している(質量保存)。
校正の結果、60cm外掛け(400 L/h)で噴流部 10–15 cm/s、平均 0.6–1 cm/s、遠端は 0.1 cm/s 未満となった。外掛けフィルターの水槽で対角の隅が滞留するという観察と整合する。投げ込み式は気泡プルームにより平均 2–3 cm/s、スポンジおよび底面式は 0.1 cm/s 以下である。
格子が粗いため数値拡散により循環は実物より弱く出ていると考えられる推定。
11数値実装上の誤り
収支が閉じるという制約は、検出手段として有効に機能した。以下、発見順に記録する。
CO2 が 3.5×1011 mg/L に発散
修正Patankar で解いた後、各プールの比 φ = cn+1/cn を用いて生物膜量や吸着量を更新している。あるプールが減衰して非正規化数(6×10−312)に達したとき、この比の計算でガウス消去が精度を失い、φ が発散した。1e-18 未満のプールを 0 に丸め(端数は台帳の除去項に計上)、比の分母に 1e-24 のガードを設けて解決した。
ソイルが窒素を生成していた(v0.7 以前から存在した誤り)
ソイルの初期 NH4 溶出を外部生成項 Pext として記述していた。すなわち窒素が無から発生していた。正しくは底床の緩効性 N(DsN)から TAN への転送である。R[TAN][DsN] に付け替え、DsN の残量でクランプして修正した。台帳を毎回検証していなければ検出できなかった。
藻類の被食を二重計上
藻類側でもエビによる摂食量を減算し、エビ側でも減算していた。algaeStep は自然な増減のみを扱い、shrimpStep が藻類の現存量から直接減算するよう構造を変更し、ステップ順を 藻類 → エビ → 水草 に固定した。
エビのデトリタス食で収支が漏れる
供給側は易分解(DfN)と難分解(DsN)の双方から見積もっているのに対し、消費は全量を DfN から減算していた。out.detF / out.detS に分離し、それぞれのプールから減算するよう修正した。最終的な残差は相対 1e-10〜1e-14 である。
孵化時に稚エビの体重が不一致
稚エビの体重を対数正規分布で個体ごとにばらつかせた一方、母体から減算する量は代表値×匹数で計算していた。実際に生成された個体の総和を減算するよう修正した。
水草の枯死で窒素が増加
loss = min(loss, B*0.5) でクランプした後に loss += B と記述しており、1.5×B が失われる(= プールに 1.5B 入る)状態だった。
溶存酸素が飽和度 350% に到達
高光量の水草水槽で溶存酸素が 28 mg/L に達した。欠落していたのは気泡としての脱気(pearling)である。葉面に付着した気泡が浮上する分を、飽和度125%以上での追加脱気項として導入したところ、150%程度で頭打ちとなった。水草水槽で日中に気泡が発生する現象に対応する。
12検証:60本のテスト
すべて Node 上で実行される(描画なし)。ブラウザでの確認は最終段階の目視のみであり、挙動の検証はすべてこちら側で行っている。
| # | 検証内容 |
|---|---|
| 1–7 | 物理化学:DO飽和(Benson–Krause)、NH3 の pKa、pH–KH–CO2 の関係、熱収支と蒸発、ヒーター故障時の冷却時定数、冷却ファンの効果 |
| 8–13 | 立ち上げ:フィッシュレスサイクリング、過剰投入による遊離アンモニア阻害、成熟槽の群集組成、濾材の選好、換水と硝酸の解析解との一致、N台帳 1e-9 |
| 14–28 | 生体:成長曲線(オスカー1年 21–29cm・2年 27–35cm)、古い水槽症候群、低酸素耐性、停電、濾材洗浄、残留塩素、塩、絶食と噛みつき |
| 29–36 | 数値と装置:修正Patankar の正値性と保存、多層濾過の直列通過、目詰まり、活性炭、ゼオライト、牡蠣殻 |
| 37–43 | 運動と流れ:遊泳運動学、個体差の CV、衝突判定(貫通ゼロ)、機材配置の妥当性、水流ソルバの発散と安定性、トレーサ粒子 |
| 44–51 | 海水:Lueker の pK1/pK2、塩分ドリフトと ATO、スキマー、ライブロック、石灰化とドージング、白化、海水の硝化 |
| 52–54 | 骨格:椎骨数と重心位置、頭部先行の転回、進行波の包絡と Videler の比較 |
| 55–60 | 水草水槽:水草の成長、日周変動、藻類の競合、エビの個体群、エビの運動学、水草水槽の N 台帳 |
テスト43(トレーサ粒子の捕捉率)のみ確率的であり、まれに不合格となる。期待値 Q·n/V に対して 0.4–2.4 倍という許容幅を設けているが、それでも統計的な揺らぎが残る。決定論的でないテストとして明示的に残している。
13未解決:シェーダの移植性
公開版には、一部の Android 端末で3D描画のシェーダがコンパイルに失敗するという未解決の不具合がある。
Assertion failed: GVI && "cannot compute gv size for oob (no global info)"
これは three.js のメッセージではなく、端末側の GPU ドライバまたは ANGLE のシェーダコンパイラが出力しているアサートである。開発環境(SwiftShader)では5槽すべてで再現しない。
GLSL ES 1.00 において必須機能でない構文に依存している箇所を6箇所まとめて修正した版を一度公開した。uniform 配列の添字をループ変数(定数添字式)のみに限定し、uniform 数を 64 ベクトルから 24 へ半減し、断片シェーダの precision highp float; 宣言を撤去し(ES 2.0 では断片シェーダの highp は任意機能であり、非対応端末ではこの宣言自体がコンパイルエラーになる)、リンク失敗を検知して簡易シェーダへ退避する経路を追加した。
結果は「魚が表示されなくなった」という報告であり、修正前より悪化したため差し戻した。
14限界
- 疾病が未実装。白点病・ネオン病・カラムナリスのいずれも扱っていない。ストレス指標 S は外見の変化と食欲低下にのみ作用する
- 水上葉/水中葉の転換がない。実際には購入直後の株がこの過程で一度崩れる
- アレロパシーがない。藻類との競合は光と栄養の取り合いのみで表現している
- 推力が運動学的。尾から生じる力で前進速度を決定しておらず、歩容モデル(Bainbridge 則)から速度を与えている。力学的に解いているのは旋回と横振れのみである
- エビの背腹方向は剛体。矢状面の屈曲は正確だが、三次元の捻れは扱っていない
- アナモックスが未実装。海水の脱窒は従来型のみである
- 室温は正弦波で与えている。実測気象データとの連携は未実装
これらの限界を含みつつ、出力される数値は現実と大きく乖離しない範囲に収まっている。60cm ハイテク水草水槽で乾重 28–39 g、ロタラの相対成長率 0.053/日(5–8 cm/週相当)、日中の DO 11.6 mg/L(飽和度150%)・夜間 8.1 mg/L、CO2 25–32 mg/L、pH 5.7–6.1、180日でトリミングにより乾重 50 g を系外へ排出、という値が得られる。
参考文献
- Arnott, S.A., Neil, D.M. & Ansell, A.D. (1998) Tail-flip mechanism and size-dependent kinematics of escape swimming in the brown shrimp Crangon crangon. J. Exp. Biol. 201, 1771–1784.
- Bainbridge, R. (1958) The speed of swimming of fish as related to size and to the frequency and amplitude of the tail beat. J. Exp. Biol. 35, 109–133.
- Burchard, H., Deleersnijder, E. & Meister, A. (2003) A high-order conservative Patankar-type discretisation for stiff systems of production–destruction equations. Appl. Numer. Math. 47, 1–30.
- Cai, Y. & Shih, H.-T. (2019); Englund, R.A. & Cai, Y. (1999) — Neocaridina / Caridina の体節比の記載
- Domenici, P. & Blake, R.W. (1997) The kinematics and performance of fish fast-start swimming. J. Exp. Biol. 200, 1165–1178.
- Kim, J.-G. et al. (2016) Hydrogen peroxide detoxification is a key mechanism for growth of ammonia-oxidizing archaea. PNAS 113, 7888–7893.
- Lighthill, M.J. (1971) Large-amplitude elongated-body theory of fish locomotion. Proc. R. Soc. Lond. B 179, 125–138.
- Lueker, T.J., Dickson, A.G. & Keeling, C.D. (2000) Ocean pCO2 calculated from DIC, TA, and the equations for K1 and K2. Mar. Chem. 70, 105–119.
- Maberly, S.C. & Madsen, T.V. (1998) Affinity for CO2 in relation to the ability of freshwater macrophytes to use HCO3−. Funct. Ecol. 12, 99–106.
- Madsen, T.V. & Cedergreen, N. (2002) Sources of nutrients to rooted submerged macrophytes growing in a nutrient-rich stream. Freshw. Biol. 47, 283–291.
- Martens-Habbena, W. et al. (2009) Ammonia oxidation kinetics determine niche separation of nitrifying Archaea and Bacteria. Nature 461, 976–979.
- Merbt, S.N. et al. (2012) Differential photoinhibition of bacterial and archaeal ammonia oxidation. FEMS Microbiol. Lett. 327, 41–46.
- Nielsen, S.L. & Sand-Jensen, K. (1989) Regulation of photosynthetic rates of submerged rooted macrophytes. Oecologia 81, 364–368.
- Oliveira, A.M. et al. (2008) — カージナルテトラの急性毒性(NH3・NO2・温度・pH の LC50/LT50)Acta Amazonica
- Sand-Jensen, K. & Gordon, D.M. (1984) Differential ability of marine and freshwater macrophytes to utilize HCO3− and CO2. Mar. Biol. 80, 247–253.
- Sauder, L.A. et al. (2024) — 海水観賞魚水槽の硝化群集(comammox 不検出、AOA/AOB の拮抗)Appl. Environ. Microbiol.
- Taylor, G.I. (1952) Analysis of the swimming of long and narrow animals. Proc. R. Soc. Lond. A 214, 158–183.
- Tropea, C. et al. (2015) — Neocaridina の胚発生と温度
- Videler, J.J. (1993) Fish Swimming. Chapman & Hall.
- Webb, P.W. (1978) Hydrodynamics: nonscombroid fish. In Fish Physiology VII.
※ 一部の文献は要点のみを二次資料で確認している。数値の責任は実装者にある。パラメータの確度 A/B/C はシミュレータの UI 内にも表示される。