GB2624036A

GB2624036A - Bespoke digital twin for chemical plant control

Info

Publication number: GB2624036A
Application number: GB2216543.5A
Authority: GB
Inventors: Caldwell Roger
Original assignee: Infrasalience Ltd
Current assignee: Infrasalience Ltd
Priority date: 2022-11-07
Filing date: 2022-11-07
Publication date: 2024-05-08
Also published as: WO2024100382A1; GB202216543D0

Abstract

A digital twin 104 of a chemical plant 100 has prises a plurality of connected blocks (fig 7) each representing part of a process carried out by the chemical plant 100. For each block, mass and energy conservation equations describe how matter and energy are conserved in the block. For each block, one or more equations describe the part of the process in the block and are formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships. The digital twin 104 is configured such that, for each block, outputs of the equations are constrained to agree with the mass and energy conservation equations.

Description

BESPOKE DIGITAL TWIN FOR CHEMICAL PLANT CONTROL

BACKGROUND

[0001]There is ongoing need to produce chemical products such as sodium bicarbonate, ammonia, metformin or others in energy efficient ways whilst still producing high quality chemical products. Process control operations may be used to drive better utilization of chemical plants and improve quality of chemical products. Process control operations typically collect data from various equipment in a large industrial plant and display the collected data at a central control room staffed by one or more human operators. This type of approach is labor intensive and subject to human error. Many factors have to be taken into account such as safety, efficiency, regulation, quality control and more.

[0002]The embodiments described below are not limited to implementations which solve any or all of the disadvantages of known chemical plant control systems.

SUMMARY

[0003]The following presents a simplified summary of the disclosure in order to provide a basic understanding to the reader. This summary is not intended to identify key features or essential features of the claimed subject matter nor is it intended to be used to limit the scope of the claimed subject matter. Its sole purpose is to present a selection of concepts disclosed herein in a simplified form as a prelude to the more detailed description that is presented later.

[0004]A digital twin of a chemical plant comprises a plurality of connected blocks, each block representing part of a process carried out by the chemical plant. For each block, mass and energy conservation equations mathematically describe how matter and energy are conserved in the block. For each block, one or more second equations describe the part of the process in the block and are formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships. The digital twin is configured such that, for each block, outputs of the second equations are constrained to agree with the mass and energy conservation equations. The digital twin may be used for controlling the chemical plant.

[0005]Many of the attendant features will be more readily appreciated as the same becomes better understood by reference to the following detailed description considered in connection with the accompanying drawings.

DESCRIPTION OF THE DRAWINGS

[0006]The present description will be better understood from the following detailed description read in light of the accompanying drawings, wherein: FIG. 1 shows a chemical plant controlled by a supervisory control and data acquisition system (SCADA) and digital twin; FIG. 2 is a flow diagram of using a digital twin to transition a chemical plant to a stable state FIG. 3 is a flow diagram of a method of using a digital twin to control a chemical plant taking into account weather measurements; FIG. 4 is a flow diagram of a method of using a response algorithm in the event of detected equipment downtime; FIG. 5 is a flow diagram of a method of working to a production target; FIG. 6 is a flow diagram of a method of designing a chemical plant; FIG. 7 is a schematic diagram of a digital twin showing blocks; FIG 8 is a flow diagram of a method of forming equations of a digital twin; FIG 9 is a schematic diagram of three blocks of a digital twin; FIG 10 is a schematic diagram of equations in one of the blocks of FIG. 9; FIG. 11 is an example of a first principles framework of a digital twin, FIGs 12 to 16 are enlarged portions of FIG. 11; FIG. 17 is an example of a phenomenological, empirical and machine learning framework corresponding to the first principles framework of FIG. 11 and for use in a digital twin such as that of FIG. 1; FIGs 18 to 22 are enlarged portions of FIG. 17; FIG. 23 shows a computing-based device in which a digital twin is implemented.

Like reference numerals are used to designate like parts in the accompanying drawings.

DETAILED DESCRIPTION

[0007]The detailed description provided below in connection with the appended drawings is intended as a description of the present examples and is not intended to represent the only forms in which the present examples are constructed or utilized. The description sets forth the functions of the examples and the sequence of operations for constructing and operating the examples. However, the same or equivalent functions and sequences may be accomplished by different examples.

[0008]The present technology uses a digital twin of a chemical plant in order to simulate the chemical plant and use the simulation results to control the chemical plant. In this way numerous efficiencies and benefits are achievable. Resources are saved such as energy or chemical inputs to the chemical plant (the chemical plant may be automated). Efficiencies are gained such as through reduced downtime of the automated chemical plant or reduced operating time of the automated chemical plant.

[0009]FIG. 1 shows a chemical plant 100 controlled by a supervisory control and data acquisition system (SCADA) 102 and digital twin 104. The chemical plant is physical infrastructure comprising reactor vessels, pipes, valves and other equipment. The chemical plant has a plurality of sensors embedded in the chemical plant which measure conditions in the chemical plant and states of equipment in the chemical plant. In some cases there are sensors in an environment of the chemical plant which measure environmental conditions, inputs to the chemical plant and outputs from the chemical plant. The sensors send measurement data outputs to the SCADA using wired or wireless connections in an automated way in real time. A non-exhaustive list of examples of sensors which are usable is: pressure sensor, temperature sensor, light sensor, acidity sensor, liquid flow rate sensor, gas flow sensor, density sensor, conductivity sensor, particle size sensor, chemical composition sensor, mass sensor, level sensor.

[0010]The chemical plant receives inputs comprising resources 108 and produces output which is a chemical product 106 and in some cases by products. The chemical plant 100 experiences weather conditions 110 and in some cases receives measurements of the weather conditions 110. The digital twin 104 receives weather forecast data 114 in some examples.

[0011]In various examples, the chemical plant is automated and is controlled by a SCADA 102. The SCADA is connected to the chemical plant 100 by wired or wireless communication and is computer implemented. The SCADA 102 comprises or is in communication with a digital twin 104 of the chemical plant 100. The SCADA and digital twin 104 are deployed in whole or in part in the cloud in some cases; thus the SCADA in some cases is a distributed control system (DCS). In other cases the SCADA and digital twin are deployed on-premises, that is, at the chemical plant [00. In some cases a programmable logic controller (PLC) is used; that is in some examples the chemical plant is controlled by one or more of a SCADA, a PLC, a DCS. [0012] The SCADA receives inputs 112 such as from a human operator via a user interface. The SCADA optionally receives other types of inputs 112 such as weather forecast data, product production target data, forecast events data. The SCADA receives measurement data from the sensors in the chemical plant and its environment. The SCADA cleans the measurement data and stores the measurement data and provides the cleaned measurement data to the digital twin. The SCADA is able to query the digital twin at any point in time to obtain values of parameters of the digital twin and other data from the digital twin. The SCADA comprises one or more algorithms, thresholds, criteria or rules which are executed to query the digital twin and use the query results to determine and send control instructions to the chemical plant 100. The SCADA is able to display query results from querying the digital twin at a user interface and to send the query results to other entities such as a mobile phone of a human operator.

[0013]The digital twin 104 is a software simulation of the chemical plant 100 and process in the chemical plant 100. The digital twin 104 has a plurality of connected blocks, each block representing part of the process carried out by the automated chemical plant 100. For each block there are mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block. For each block there are one or more second equations describing the part of the process in the block. The second equations are formed using one or more of first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships. The digital twin is configured such that, for each block, outputs of the second equations are constrained to agree with the mass and energy conservation equations at a given time step. This gives a highly accurate digital twin which produces simulation results for controlling the chemical plant. The digital twin gives a more accurate simulation than has previously been possible in practical time scales. The accuracy is achieved in this practical manner by using a large number of second equations (since there are a plurality of second equations for each block) and by constraining the second equations to agree with the mass and energy conservation equations at each time step.

[001 4]First principle relationships are equations or terms in equations which express the laws of physics, thermodynamics, chemistry or other laws of the physical universe. A non-limiting example of a first principle relationship is a stoichiometric chemical reaction described using a chemical equation that expresses the elementally balanced conversion of one set of molecules to another set, including the distinction between subspecies and phase changes where appropriate. [001 5]Phenomenological relationships are explained in more detail below as are machine learning derived relationships.

[001 6]The digital twin operates in a first time interval such that the outputs of the second equations are constrained to agree with the mass and energy conservation equations in the first time interval. The method is repeated for subsequent time intervals by computing the outputs of the second equations at each time interval. In this way the digital twin simulates ongoing operation of the chemical plant 100 in an accurate manner because, at each time interval, the second equations are constrained to agree with the mass and energy conservation equations. [001 7]As mentioned above, the second equations are formed using phenomenological relationships in some cases. Phenomenological relationships are relationships expressed mathematically that take into account physical attributes of the specific apparatus of the chemical plant 100 and/or chemical attributes of the process in the chemical plant 100. By taking into account physical attributes of the chemical plant 100 and/or chemical attributes of the process in the chemical plant 100 it is possible to simulate the chemical plant and its process accurately. A phenomenological relationship is one which may change based on the characteristics of the chemical plant such as size of materials whereas in contrast first principles relationships are invariant to characteristics of the chemical plant. By having the equations with phenomenological relationships constrained to agree with the mass and energy conservation equations it is possible to find accurate solutions in practical time scales since the number of possible solutions is reduced. A non-exhaustive list of examples of physical attributes of a chemical plant which are expressed mathematically using phenomenological relationships is: dimensions of vessels of containment, material characteristics of vessels of containment (such as coefficients of friction, thermal conductivity, elasticity, sound reflectivity, light transparency, external vibrational forces, radiation absorption, porosity, permeability, electrical conductivity), internal structural elements and their respective characteristics, energy sinks and transmitters On the forms of heat, electricity, electromagnetic radiation), gas diffusers, liquid injectors, spray nozzles, coatings, and catalytic surfaces. A non-exhaustive list of examples of chemical attributes which are expressed mathematically using phenomenological relationships is: atoms or molecule concentrations, heat, catalytic surfaces, ionic strength, radiation absorption, porosity, permeability, electrical conductivity, electromotive forces, electricity, electromagnetic radiation. Such chemical attributes are used in phenomenological relationships in the digital twin to describe chemical bonding of atoms and molecules (such as atoms, covalent, ionic, neutral, anionic, cationic species) to form new molecules.

[0018]The second equations are formed, in some cases, from first principle relationships and applied (to the simulation of a particular chemical plant) by lumping equation parameters together as tunable constants. This gives a principled way to obtain the second equations which facilitates accuracy. The empirical relationships mathematically relate inputs and outputs of a rate process in a block of the digital twin and are developed and adjusted periodically using error regression comparison of predictions from the digital twin 104 and corresponding measurements from the chemical plant 100. In this way accuracy is achieved since the empirical relationships describe observed rate processes. By having periodic adjustment it is possible to maintain accuracy despite drift in sensor measurements, wear and tear of equipment and changes in the rate processes.

[0019] As mentioned above, the second equations comprise machine learning derived relationships in some cases. In this case values of lumped equation parameters are predicted from appropriately trained machine learning systems. Any suitable conventional machine learning system is used such as a neural network which takes as input measurements from the chemical plant, and/or parameter values derived from earlier blocks of the digital twin, and produces as output a prediction of a value of the lumped equation parameters. The machine learning system is trained using supervised machine learning in a preferred example since this gives high accuracy and since it is possible to obtain large amounts of labelled training data. The labelled training data is obtained by recording measurements from the chemical plant 100 and calculating corresponding ground truth values of the lumped equation parameters using first principles relationships. Note that it is not essential to use supervised machine learning and that in some cases, such as where it is difficult to compute ground truth values of the lumped equation parameters, unsupervised or semi-supervised learning is used.

[0020]In some cases the machine learning derived relationships are adjusted periodically by retraining the machine learning system using recently obtained training data. By having periodic adjustment it is possible to maintain accuracy despite drift in sensor measurements, wear and tear of equipment and changes in the rate processes [0021]The plurality of blocks comprise: at least one block that has a representation of a thermal boundary to the ambient environment and at least one block that comprises a representation of one or more of the following processes: a stoichiometric chemical reaction, a dissolution from solid to liquid phase, a precipitation from liquid to solid phase, an absorption from gas to liquid phase, a vaporization from liquid to gas phase, separation of material by phase of matter, separation by centripetal forces, separation by filtration.

[0022]In this way the digital twin is versatile as it is able to represent a variety of types of process. By connecting a plurality of blocks which represent different processes the versatility is further enhanced. The modular nature of the blocks facilitates their reuse between different digital twins. Having at least one block that has a representation of a thermal boundary to the ambient environment is found to be particularly important for improving accuracy of control of the chemical plant. It is also found that taking into account weather forecast data 114 in the digital twin 104 gives significant performance improvements at the chemical plant 100.

[002 3]In various examples, the process carried out by the automated chemical plant is any of: production of sodium bicarbonate, production of ammonia, production of metformin, desulfurization of flue gas, carbon dioxide removal from a gas using an amine solvent, the production of hydrogen using solar power. Using the digital twin it is possible to achieve higher efficiencies in these chemical processes.

[0024]FIG, 2 is a flow diagram of using a digital twin 104 to transition a chemical plant 100 to a stable state. An output of the simulation in the digital twin 104 is used to change at least one controllable input to the chemical plant 100 in order to compensate for uncontrolled variations in at least one input to the chemical plant 100. In an example, the input to the chemical plant 100 is a stream of gas where a carbon dioxide CO2 concentration in the stream is measured 200. The measured CO2 concentration values are sent to the SCADA 102 which updates 202 the digital twin 104. The CO2 concentration in the input gas varies due to uncontrollable conditions in gas supply.

[0025] The SCADA 102 receives predicted parameter values from the updated digital twin 104 (such as a predicted quality of the product 106) and uses these to determine whether to maintain quality of the product 106 or not at decision point 204. The decision at decision point 204 is computed by the SCADA 102 using one or more rules, thresholds or criteria.

[0026]In response to a decision to maintain quality at decision point 204, the total gas flow is changed 212 to keep the overall reaction rate constant. When the digital twin 104 predicts that this change is not enough the decision at decision point 204 is to not maintain quality and to transition 206 the chemical plant to a stable state. In this case the SCADA controls the chemical plant in order to keep the reaction stable for a potentially deteriorating but foreseen condition. In this condition, for example in a sodium bicarbonate production plant where the CO2 availability has dropped, a reduction in carbonate flow rate in combination with an increase in total gas flow is made to keep output parameters, such as super saturation, from changing rapidly. In another example where the chemical plant 100 is producing pure CO2 106 from a flue gas stream using an amine solvent, the response may be to increase the flue gas flow rate to keep the production of CO2 constant.

[0027]In the event that the uncontrolled change continues to deteriorate, a check at check point 208 determines there is a predicted emergency condition. The SCADA implements a set of emergency actions to ensure that the chemical plant 100 moves to a safe state 210 of operation (such as partial shutdown or complete shutdown).

[0028]In response to the simulation predicting that a quality of output of the chemical plant can be maintained, the SCADA changes the at least one controllable input; and in response to the simulation predicting that a quality of output of the chemical plant cannot be maintained, the SCADA transitions the chemical plant to a stable state; and in response to the simulation predicting continued deterioration of the uncontrolled variations, the SCADA moves the chemical plant to a stable state.

[0029]FIG. 3 is a flow diagram of a method of using a digital twin to control a chemical plant taking into account weather measurements. Weather measurements are received 300 from sensors in or around the chemical plant 100. In some cases the weather measurements are replaced or supplemented with predicted weather measurements from a weather forecasting service in communication with the SCADA. The weather measurements and/or predicted weather measurements are used to update 302 the digital twin and the updated digital twin is then used by the SCADA to control 304 the chemical plant 100. Here the digital twin 100 has at least one block that has a representation of a thermal boundary to the ambient environment. Because the ambient environment (i.e. the environment outside and around the chemical plant rather than within the chemical plant) is taken into account there is significant improved accuracy of the digital twin and thus also of the chemical plant control. For example, the simulation in the digital twin 104 is used to predict variations in reactions due to weather conditions and the predictions are used to control the chemical plant 100 to compensate for excursions in the process that would otherwise be caused by the changes in the weather. [0030]In an example the digital twin comprises equations describing heat transfer through walls of pipes and vessels in the chemical plant, and the digital twin is configured to predict the temperature of a primary reaction zone in the chemical plant from weather conditions as well as from measurements obtained from a physical temperature probe in the chemical plant but outside the primary reaction zone. This is extremely useful where it is not possible to have a physical temperature probe (i.e. temperature sensor) in the primary reaction zone due to extremely harsh physical and chemical conditions in the primary reaction zone. By using the digital twin to predict the temperature of the primary reaction zone from weather data in combination with the closest available temperature probe data significant improvements in accuracy are obtained. [0031]In some cases the predicted temperature of the primary reaction zone is used to control an external heat sink on an input stream to the chemical plant.

[0032]Phenomenological equations are used to describe heat transfer through the walls of pipes and vessels, and ambient weather forecasts are used to describe the effects of the environment. Mitigation methods include adjustments for critical temperatures whose accuracy cannot be measured directly. For example, the temperature probe for a primary reaction may be located some distance away (such as near the wall of a vessel) from the primary reaction zone (such as in the center of the reactor and not reachable with an instrument). As part of the mitigation calculations, the digital twin calculates the theoretical temperature of the primary reaction zone and compares it to the nearest temperature probe. It has been discovered that the difference between the actual temperature and the nearest temperature probe correlate to changing ambient temperature differences, and that these differences cannot be regressed; hence; a mitigation algorithm based on changing weather is introduced to compensate.

[0033]One form of a weather based mitigation algorithm is to take the weather forecast high and forecast low daily temperatures and introduce a forced change to the measured reactor temperature, which is implemented through an external heat sink on one of the input streams (for example, through a slurry recycle stream) to compensate for the difference between the maximum and minimum differences between the measured reactor temperature and the simulated primary reactor zone temperature.

[0034]If weather forecast data is available, this mitigation algorithm is adjusted to include more detail in the ambient forecasts, including hourly, windspeed, precipitation, humidity, cloud cover, and radiation forecasts. In this way significant performance improvements are found In some examples, the mechanical design of the chemical plant includes a heating jacket on a reactor vessel. Using a heating jacket gives more control over the environmental effects, especially in more extreme climates, as the digital twin is able to make setpoint changes to the jacket temperature, in addition to a heat sink.

[0035]FIG. 4 is a flow diagram of a method of using a response algorithm in the event of detected equipment downtime. The SCADA receives status data from equipment in the chemical plant 100 and is able to detect 400 equipment downtime. Equipment downtime occurs where equipment fails or is scheduled for maintenance. Inputs to the SCADA from a human operator include scheduled maintenance of equipment in some cases.

[0036]The SCADA detects 400 downtime of equipment in the chemical plant, selects 402 a response algorithm from a store 404 of response algorithms according to characteristics of the detected downtime of equipment, and applies 406 the selected response algorithm to control the chemical plant in combination 408 with predictions from the simulation.

[0037]On occasion ancillary equipment can fail (or be scheduled for maintenance) that has an effect on the process in the chemical plant. In this situation, the simulation of the digital twin is used to mitigate the downtime by having a preset response algorithm depending on the types of downtime situation. The store of response algorithms is developed in advance for any of the conceivable equipment failures or scheduled maintenance procedures. Using the store of response algorithms significantly reduces human error, making the chemical plant 100 more efficient and safe to operate. Using the digital twin 104 to control the chemical plant in parallel with the selected response algorithm is found to give significantly improved performance as compared with using the response algorithm alone.

[0038]FIG. 5 is a flow diagram of a method of working to a production target.

[0039]The SCADA receives 500 a target for production of an output (the chemical product and/or by product) by the chemical plant 100 and, based on a difference between a predicted output computed by the simulation with the target, controlling the chemical plant 504 so as to reduce the difference.

[0040]In an example, it is desired for the chemical plant 100 to achieve production targets which vary from day to day based on orders and inventory. Weather forecasts and other known variable inputs are programmed to give a forecast output of the chemical plant 100 based on preset controlled variable inputs, such as input gas flow and carbonate feed flow. The simulation results from the digital twin are used by the SCADA to control the set points of the inputs 108 to the chemical plant to achieve variable overall reaction rates (production rates) to achieve production targets with very low variation. This allows the chemical plant to run at the speed required to meet the production target and avoids building excessive inventory.

[0041]The target in operation 500 of FIG. 5 is any one or more of a quantity of the output, a purity grade of the output, a grade of particle size distribution of the output.

[0042]Product grade distribution demand may shift (for example grades of product purity, or grades of particle size distribution) from time to time based on business conditions. The simulation can be used to manipulate input conditions to achieve the grades desired, with known impact on production times determined in advance.

[0043]Reactor residence time can affect particle size, and the simulator predicts the changes needed (such as feed concentration) to atTect the final particle size distribution needed to achieve the customer order specifications.

[0044]Higher purity products may require running lower overall reaction rates to achieve higher purity and tighter particle size distribution. These impacts on productivity are calculated by the simulator in advance and become a critical tool for optimizing plant productivity.

[0045]The digital twin and SCADA are used for early detection of faulty equipment such as miscalibrated instruments or changes in pump performance. This is done by dynamically comparing a simulated result with a measured result.

[0046]In an example, flow rate calibration is checked by comparing a measured flow rate to a simulated flow rate from the digital twin. Small differences show up as an accumulation in a simulated volume resulting from the flow rate, indicating a measurement error in a flow meter. Even if the errors in the suspected meters are small, the errors can be tracked and equipment scheduled for calibration/maintenance prior to a process failure. This type of mathematical analysis is extremely sensitive and the calibration errors are almost impossible to detect by any other known method.

[0047]In another implementation the SCADA is programmed to dynamically compare a rate measured in the chemical plant with a simulated rate from the digital twin, and in response to the comparison meeting criteria, detects an error in a flow meter in the chemical plant.

[0048] In another implementation example using the digital twin the SCADA is programmed to compare a measured temperature of a process stream in the chemical plant after a heat exchanger to a temperature simulated by the digital twin; and detects scaling on surfaces of the heat exchanger based on the comparison.

[0049] The digital twin is used to carry out regression of plant wear and tear. It is normal for parts and instruments in the chemical plant 100 to have deviations over time that create error and noise in output signals needed to maintain control. As described above, empirical relationships are used in the second equations within the digital twin. The empirical relationships comprise empirical equations developed using error regression that result in a set of constants which the digital twin uses that minimize error and improve predictive capability based on mathematical statistics. As the plant wears over time, the regression is repeated using data from the more recent time frame, thus recalibrating the digital twin to the current conditions, improving its predictive capability.

[0050]By comparing the updated regression correction factors to the predecessor correction factors, changes that indicate progressive deterioration can be identified giving insight to equipment that may be failing faster than others. This information is forwarded to the SCADA which can be programmed to monitor changes in error regression correction factors in order to detect equipment failing faster than other equipment in the chemical plant.

[0051]FIG. 6 is a flow diagram of a method of designing a chemical plant. A design is created manually or using an automated process and physical characteristics 600 of equipment according to the design are determined. The physical characteristics are input to the digital twin in order to set parameter values of one or more of the equations in the digital twin. The digital twin is then executed 602 and the SCADA queries values of parameters in the digital twin and compares them with one or more thresholds at check point 604. If the thresholds are met the design is deployed 606. If the thresholds are not met the design is adjusted 608 and the process repeats. [0052]Thus the digital twin is usable to predict performance of any one or more of: reactor vessels of different dimensions, gas sparging apparatus, insulation, physical characteristics of the chemical plant. The digital twin is used to predict performance of the chemical plant for a range of equipment performance parameters of equipment in the chemical plant. The digital twin is also used to predict performance of the chemical plant for a range of weather conditions. [0053]The simulation can be used to find optimal vessel dimensions, and other parameters, before the vessels are built. Reactor vessel diameter and height is simulated in order to ensure good performance over a range of operating conditions that are predicted. The same is true for other components and dimensions, such as gas sparging apparatus, insulation, and a wide variety of phenomenological parameters. By running the simulation across the range of anticipated process conditions, the optimum nominal physical dimensions of equipment can be defined. [0054]Once specific equipment is chosen (such as specific pumps and other off the shelf capital), the simulation is used to ensure process performance across the range of specified equipment performance parameters. The output of a chemical plant is simulated across a full range of anticipated conditions, such as annual weather forecast, to ensure that the design will remain efficient under the conditions of a specific site.

[0055]FIG. 7 is a schematic diagram of a digital twin 104 showing blocks. In the example of FIG. 7 there are three blocks connected in series although it is possible to use more than three blocks or only two blocks. As mentioned above, each block represents part of the process being simulated by the digital twin 104 and each block comprises a plurality of mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block. For each block there are one or more second equations describing the part of the process in the block. The second equations are formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships. The digital twin is configured such that, for each block, outputs of the second equations are constrained to agree with the mass and energy conservation equations at a given time step. Therefore using blocks gives the benefit of improved efficiency and accuracy of the digital twin.

The blocks are modular by nature and are reusable between different digital twins. The choice of the number of blocks to use in a particular digital twin 104 is made depending on a variety of factors such as the architecture and topology of equipment in the chemical plant and to cover the key measurement and control points covered by the digital twin.

[0056]FIG. 8 is a flow diagram of a method of forming equations of a digital twin. The process of the chemical plant is first divided into parts, one per block. The process of FIG. 8 is then followed separately for each block.

[0057]Given a block as input, the automated method begins by defining 800 first principles equations of the block with the overall reaction rate embedded therein. The definition for the overall reaction may be received as input from an operator. The overall reaction is then broken into substeps by using rules or other criteria to determine the substeps. For each substep, a rate equation is specified 804 by selecting an appropriate rate equation from a database. A framework is then stored comprising a plurality of rate equations 806.

[0058]The second equations are then formed. The second equations are formed using one or more of first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships. As shown in FIG. 8 phenomenological equations are formed 808 by mathematically describing physical attributes 812 of the apparatus in the chemical plant, and/or the physical and chemical attributes of molecular behaviour, and by lumping equations parameters together as tunable constants. A check is made at decision point 810 as to whether the phenomenological equations are sufficient in terms of various criteria. The check is automated and is made using rules or thresholds. The criteria are one or more of: accuracy, computational efficiency, amount of available measured data, fidelity of measured data, size of dataset needed for machine learning.

[0059]If the phenomenological equation is determined to be sufficient for use in the digital twin it is stored 814.

[0060]if the phenomenological equation is determined to be insufficient for use in the digital then either an empirical equation 816 or a machine learning model 818 is used in place of the phenomenological equation.

[0061]In cases where the understanding of the particular physical or chemical rate phenomenon is insufficient or where there is not enough measured data available or where the measured data is insufficient in its fidelity or requires first principles descriptions at dimensional scales where computation effort becomes impractical, then equations expressing empirical relationships are used to mathematically relate the required inputs and outputs of these rate processes. The empirical rate equations are algebraic, differential or partial-differential equations.

[0062]An alternative approach is to utilize machine learning and data science disciplines to recognize and describe these rate processes. Machine learning models are used where there is ability to collect large amounts of training data.

[0063]The machine learning models and empirical equations are stored for use by the digital twin. The machine learning models and empirical equations are optionally updated by retraining the machine learning models using recently obtaining training data, or by regressing the empirical equations using recently observed measurements. The updates are done at regular intervals and/or when performance of the digital twin is found to be below a threshold. [0064]FIG. 9 is a schematic diagram of three blocks of a digital twin in the case that the chemical plant 100 is producing sodium bicarbonate. A first one of the blocks (Block-1 in FIG. 9) represents a sodium bicarbonate reaction chamber which carries out the following chemical reaction: Na2CO3(a) + CO2(g) + H20(1) 4 2NaHCO3(s) [0065]which is expressed in words as, sodium carbonate in aqueous form, plus carbon dioxide in gaseous form, plus water in liquid form produces sodium bicarbonate in solid form. [0066]In the example of FIG. 9 the digital twin comprises three connected blocks being a first block representing a reaction chamber process, a second block representing a slurry recycle, heat sink and production offtake process, and a third block representing a gas recycle and makeup carbon dioxide (CO2) gas process.

[0067]Block-1 receives an aqueous feed input which is a liquid feed input consisting of an aqueous solution containing sodium carbonate (SC), sodium bicarbonate (SBC), and water (H20), whereby the concentrations, temperature, and flow rate are measured with instruments in the chemical plant 100 that are common in the art and sent to the SCADA.

[0068]Block-1 also receives a gas input boundary from Block-3 containing CO2, H20, and other gases that are inert, whereby the behavior of the gases are modeled to be ideal. The concentration of CO2, temperature, total mass flow rate, and pressure of the input gas are measured with instruments in the chemical plant that are common in the art and sent to the SCADA.

[0069]Block-1 represents a cylindrical vessel (reactor) with a specified diameter (D), and height (H), whereby the input gas enters at the bottom through a set of orifices of small diameter. The gas is modeled as a mixture of bubbles with a mean radius (R), whereby the total gas / liquid interface area (A) is described as an area equivalent to the volume of gas divided by the volume of an average bubble size times the surface area of a bubble of average size. Measurement of volume of the condensed (solid-liquid-bubble) phase in the vessel is measured by physical sensors and sent to the SCADA. Measurement of pressure at the base of the cylinder and at the top of the cylinder is measured with instruments that are common in the art and sent to the SCADA.

[0070]Block-1 also has an output boundary to Block-2, wherein a slurry mixture of solid and liquid mass is output from Block-2, said slurry containing aqueous SC, aqueous SBC, and solid SBC crystals, wherein the temperature, density, pH, conductivity, aqueous SC concentration, aqueous SBC concentration, total mass flow rate, and solid crystal size distribution are measured with instruments that are common in the art and sent to the SCADA.

[0071]Block-1 also has an input boundary from Block-2, wherein a slurry mixture is received containing aqueous SC, aqueous SBC, and solid SBC crystals, wherein the temperature, density, pH, conductivity, aqueous SC concentration, aqueous SBC concentration, total mass flow rate, and solid crystal size distribution are measured with instruments that are common in the art and sent to the SCADA.

[0072]Block-1 has a thermal boundary to the environment, wherein heat is transferred from Block-1 to the environment outside the vessel through a surface determined by the outer surface area of the vessel.

[0073]Block-1 has an output boundary to Block-3, wherein a gas mixture containing CO2, H20, and other gases that are inert is output to Block-3. The behavior of the gases is modeled to be ideal. The concentration of CO2, temperature, flow rate, and pressure of the output gas is measured with instruments that are common in the art and sent to the SCADA.

[0074]Block-2 represents slurry Recycle, heat sink, and production otliake.

[0075]Block-2 has a mixed phase input received from Block-l's output boundary as described above. Block-2's solid output is removed as a product offtake, wherein characteristics are identical to mixed phase input from Block-1, except that mass flow rate is a fraction between 01.0 of the flow rate from Block-1.

[0076]Block-2's mixed phase output is removed through the boundary with Block-1 and where the mass flow rate to Block-1 is the difference between the inflow rate from Block-2 to Block-1 and Block-2's product offtake flow rate.

[0077]Block-2 represents a physical apparatus in the chemical plant such as stainless steel piping of specified diameter, and length, wherein portions of the outer surface of the return line to Block-I is exposed to an annulus of cooling water. The total flow rate through the piping is driven by a slurry pump, and the fraction of slurry that is diverted to product offtake is metered with a control valve.

[0078]Block-2 has a thermal boundary to cooling water, wherein heat is transferred to a cooling water stream, wherein the flow rate, input temperature, and output temperature of the cooling water and the slurry stream is measured with instruments that are common in the art and sent to the SCADA.

[0079]Block-2 has a thermal boundary to the environment, wherein heat is transferred from Block-2 to the environment through a surface determined by the outer surface area of the piping that is not exposed to the cooling water.

[0080]Block-3 represents a gas recycle and makeup CO2 gas process.

[0081]Makeup CO2 is provided from an external source to Block-3, wherein the temperature, pressure, and mass flow rate of CO2 gas is measured with instruments that are common in the art. Makeup CO2 is carbon dioxide introduced in order to maintain pressure of gas in the chemical plant above a specified threshold.

[0082]Block-3 represents a physical apparatus in the chemical plant which is described as stainless steel piping, wherein gas is driven through the piping with a compressor capable of reaching the gas inlet pressure of Block-1 (high pressure), wherein the inlet to the compressor (low pressure) is supplied by the combination of makeup CO2 gas and the outlet gas from Block-1. Post compression, the gas is passed through a cooling water heat exchanger, prior to exiting Block-3 as the gas input boundary to Block-1.

[0083]Makeup gas, recycled gas (input from Block-1), and combined outlet gas (output to Block-1), all have temperature, pressure, and mass flow rate measured with instruments that are common in the art and the measurements are sent to the SCADA.

[0084]More detail about the equation development for the second equations (i.e. the first principles equations with phenomenological, empirical or data science relationships) in the digital twin is now given for the case where the digital twin simulates Block-1 of FIG. 9 or any other reaction chamber for producing sodium bicarbonate according to the equation: Na2CO3(a) + CO2(g) + H20(1) 4 2NaHCO3(s) [0085]ln this example one of the blocks of the digital twin is broken down into a plurality of sub steps. The first block (Block-1) comprises three substeps that enable discrete speciation of chemical components, each substep represented using an equation, whereby the three equations are solved simultaneously in a given time interval using the digital twin.

[0086]FIG. 10 is a schematic diagram of reactions in one of the blocks of FIG. 9 (Block-1) and shows these separated into three substeps.

[0087]The first step in developing equations is to break down the overall reaction equation into substeps that enable discrete speciation of chemical components so that physical properties of the species can be used to describe the reaction rates and thermodynamics The main SBC formation reaction in Block-1 is broken down into 3 substeps each being an equation. All of these steps happen simultaneously within Block-1; however, they must happen in sequence for a given chemical species. For example, a carbon atom in the incoming CO2 gas goes through equation 1, then equation 2, and then equation 3 in order to become a carbon atom in the sodium bicarbonate solid output. These steps are illustrated as an example in FIG. 10.

[0088]Equations are developed that describe the reaction rate of each substep, and then these equations are solved simultaneously over each arbitrary unit of time. The equations involve complex relationships between properties of the chemical species, the physical environment that they are contained within, and mathematically derived constants that are used in place of characteristics that are unknown. However, by breaking the stoichiometric chemical equations into substeps, it begins to limit the degrees of freedom for the number of solutions possible for each set of first-principles equations within a block.

[0089]For instance, a physical characteristic such as a gas bubble size in the first gas absorption step can logically be assumed to not affect the precipitation rate of sodium bicarbonate NaHCO3 from liquid to solid in the third step. Likewise, a physical characteristic such as the total crystal solid surface area can have a great effect on the precipitation rate in the third step but can logically be assumed to have almost no effect on the first gas absorption step. By subdividing the block into physical transformations which each are uniquely described, the degrees of freedom for the simultaneous solution for all of the equations is reduced. By increasing the number of phenomenological equations embedded in the first-principles mass and energy conservation equations, the degrees of freedom continue to be constrained which reduces the reliance on mathematically derived constants to solve the simultaneous equations, therefore increasing the net confidence in the simulation's predictive results.

[0090]Fundamentally, the finer the resolution of the first principle equations that define a block, the greater the confidence in the predictive results of the simulation. In a block containing a stoichiometric chemical reaction, this means that the more precise the understanding of the sequence of sub-steps in an overall reaction, the more reliable the predictive capability of the simulation will be for the very simple reason that the universe of optional solutions is decreased. But there are very significant limitations to the resolution that can be obtained because of the inability and inaccuracy of the instruments that are used to measure the dynamic conditions. Therefore, there are always trade-offs in this method of hybrid algorithms that contain both fundamental and empirically derived terms. One limit is determined by computational capacity, and another limit is determined by measurement capability.

[0091]In contrast to previous approaches the digital twin comprises a large array of fundamental and empirical terms used in a multitude of simultaneous equations to find simultaneous solutions over an arbitrarily short duration of time. This approach is found to give improved accuracy. [0092]More detail about how the digital twin simulates the first substep of FIG. 10 is now given. The first substep represents absorption/desorption of gas/liquid (Step-1 physical step). The first substep represents the physical absorption of carbon dioxide gas into an aqueous liquid, which is performed in Block-1 [0093]As explained with reference to FIG. 8, equation development involves defining mass and energy conservation equations for a block. In the present example this is done for a substep. The rate of CO2 absorption in a reaction chamber block, in this case a bubble column, is dependent on the physical attributes of the gas dispersion device, for example, a rigid sparger system of gas pipes with small perforations or holes at the bottom of the column, as well as the physical attributes of the bubble column itself, for example, the column diameter and height, and the velocity components imparted by the flow rate and direction of the bulk fluid through the reaction chamber block, as well as the local fluid velocity components imparted by physical geometric arrangements, such as annulus subdividers, baffles, tangentially positioned recirculation pipes, etc., as well as the physical attributes of the fluid in the block, for example, the composition of the aqueous (aq) phase, and the mass and size distribution of solid SBC particles in the fluid mixture. The physical CO2 gas absorption rate in this three-phase (gas, liquid, particle) fluid is therefore a complex function of many variables that are also not necessarily constant over the complete block control volume. The digital twin comprises equations that describe the block mass, and energy conservation. In some cases, the mass conservation equations comprise particle population conservation equations. Particle population is expressed along a particle size coordinate. Therefore using particle population conservation equations in the digital twin allows calculation of changes in the particle size distribution as the block states change over time and adds fidelity to the digital twin. Any suitable known equations which describe block mass and energy conservation are used. Such equations are equally applicable to any sub-dimensional block within the reaction chamber block so that the mass, energy and particle population are simultaneously satisfied for each of the sub-dimensional blocks. A sub-dimensional block is part of a block where a block is partitioned into specific subsections to add more fidelity. For example, to subdivide the block into different reaction regimes, such as a control volume above a sparser where the temperature tends to be higher than in the bulk slurry.

[0094]As explained with reference to FIG. 8, phenomenological equations are formed 808 and this is done for the substep. The relationships specific to the physical absorption of CO2 in the bubble column reaction chamber block example are described in the form of phenomenological equations, where the functionalities are rooted in first principles theory but then applied by lumping equation parameters together as tunable constants. In cases where the understanding of the particular physical or chemical rate phenomenon is insufficient or where there is not enough measured data available or where the measured data is insufficient in its fidelity or requires first principles descriptions at dimensional scales where computation effort becomes impractical, then empirical relationships are used to mathematically relate the required inputs and outputs of these rate processes. An alternative approach is to utilize machine learning and data science disciplines to recognize and describe these and other sub-step rate processes.

[0095]Different examples of phenomenological and empirical equations for the physical absorption of CO2 (and the other example Block-1 substeps) are shown below, but importantly, all these equations are embedded within the above time-based (dynamic) mass, energy and particle population conservation framework that mathematically describes how matter and energy interact within arid between each process block arid substep.

[0096]The physical absorption and evaporation of molecules are expressed by their stoichiometric reactions. For example, for the carbon dioxide (CO2) molecule, absorption and evaporation is represented as follows (with analogous expressions for 02, S02, etc.): CO2(g) <=> CO2(a) (RI) [0097]Which is expressed in words as, rate R1 denotes a rate of change of carbon dioxide molecules between a gaseous state CO2(g) and an aqueous state where carbon dioxide is dissolved in water CO2(a).

[0098]The interfacial gas-liquid contact areas are i) at the slurry surface of the column and ii) between the individual bubbles and the surrounding liquid as they rise up in the column. The transfer rate across the interface layer is treated phenomenologically by using a theoretical concept, such as a film model that describes a stagnant film of small thickness at the surface of the liquid adjacent to each gas bubble. A film model follows from conducting a first-principles mass balance over the interface between a dispersed gas bubble and the surrounding liquid phase. The difference between the solubility at the gas-side of the bubble interface (CO2 eq) and the dissolved concentration at the liquid side (CO2(a)) drives the molar diffusion rate (dn/dt) across this interfacial film layer and is denoted as: dn(CO2)/dt kLa*([CO2 eq] -[CO2(a)]) (El) [0099]Which is expressed in words as, the molar diffusion rate of carbon dioxide is equal to the average interfacial mass-transfer coefficient (denoted kLa) of the bubbles multiplied by the difference between the solubility of carbon dioxide at the gas-side of the bubble interface and the dissolved concentration of carbon dioxide at the liquid side. Equation El is an example of a phenomenological equation because El has its roots in the first-principles mass conservation law around a gas bubble, but has been formed by making a pseudo steady-state assumption, followed by lumping the average linear diffusion coefficient through the stagnant interfacial film layer (1th in units of m/s metres per second) and the average total bubble area per unit volume (a in units of m2/m3 metres squared or metres cubed) together to form a single lumped parameter kLa (in units of 1/s one per second), which is referred to the average interfacial mass-transfer coefficient.

[00100] The thermodynamic states of the gas and liquid phases fix the equilibrium value [CO2 eq] at time t, which is implemented either via embedded first principles algorithms, or, embedded data science algorithms developed using expert thermodynamic software in combination with, for example, machine learning tools, or, any combination thereof. Similar approaches are followed for all gaseous molecules, including for the solvent phase. For example, for the water (H20) molecule, absorption and evaporation is represented as follows H20(1) <=> 1-120(g) (R2) [00101] where the thermodynamic states fix its equilibrium vapor pressure pH20 eq so that the rate of evaporation across the interface is: dn(H20)/dt = kLa.(pH20 eq -p1120) (E2) [00102] Which is expressed in words as, the evaporation rate of water is equal to the total interfacial area multiplied by the difference between the equilibrium vapor pressure of water at the gas-side of the interface and the bulk partial pressure state of gaseous water in the specific sub-block.

[00103] Since the states at each block, or sub-block, boundary node are described by the time-based first principles mass and energy balance equations for each species, Equations El and E2, and similar equations for the other gaseous molecules, are only concerned with describing the rate of molecular transfer between the liquid and gas phases in that specific block or sub-block. Furthermore, with the equilibrium states known, as described above, and again using the above bubble column as an example, the mathematical connection with the non-equilibrium side of each interface is fully established when the physical or chemical consumption rate of that species has been calculated. For example, in the case of CO2, the value of CO2 aq is known at each time step after its chemical conversion rate (see next subsection below) has been added to the column block algorithms. Similarly, in the case of H20, the time-based value of pH20 is known at each time step when the pressure/volume state of the free board gas space above the slurry inventory volume are mimicked and, where necessary, any vapor condensation rate algorithms have been implemented, for example, at specific block boundaries, such as at the top of the column, exposed to the ambient environment outside the column block.

[00104] With the above mathematical constraints in place, the mass transfer rate across the column surface at the slurry level interface (i) is dealt with by using a lumped parameter approach to define the average interfacial mass-transfer coefficient kLa, since the interfacial area (a) is known for most block geometries. Finding a realistic value for the surface Itha is therefore possible but relatively inconsequential since the most important perturbations are related to the transients around the gas holdup and the related (ii) total bubble surface area.

[00105] The gas holdup and total bubble area in the column block is a complex function of the gas flow, physical properties of the gas dispersion device (e.g., the rigid sparger design), slurry phase properties and the column hydrostatic head. This parameter space is mapped by performing displacement, optical and other measurements in test columns, even at different hydrostatic heads. Where such information is unavailable, the digital twin is configured to predict the behaviour of the rising bubbles as functions of the key parameters. Either way, the differential mass balance equations that describe the transient behaviour of chemical species mass in the gas-liquid-solids bulk phase in the column block or sub-block, is partitioned between the condensed liquid-solids slurry phase species and the bubble gas phase, which includes an arbitrary small mass of liquid representing thin liquid films around the bubbles that would tend to move with the bubble swarm. The same interfacial mass transfer rate equations as the El and E2 examples above (or analogous equations for any volatile molecule) are used to describe the different species mass exchange rates between these sub-partitions. As with the gas/liquid exchanges at the interface between freeboard gas space and the column slurry surface (described above), the rates and directions of the interfacial mass exchanges are determined by the concentration gradients between the actual species states and their thermodynamic states at any point in time. The gas hold-up is therefore an important element of the mass balance and is set up so as to minimize the difference (error) between that calculated hold-up volume and a target value. The target value refers to empirical knowledge of how much gas volume (typically expressed as the volume fraction of dispersed bubbles) is held up at any point in time in columns of the reaction chamber as a function of the specific column operating states, for example, the gas superficial velocity, sparger type, column level, slurry density, etc. This target value is obtained by one or more of describing by fluid-mechanics algorithms (of dimensionless numbers), obtained from the literature, obtained from empirical equations or lookup tables mapped by physical measurement over ranges of typical operating points for a specific column. In some cases the target value is tunable by an operator, but usually is pre-determined for specific equipment via dedicated testing. The gas hold-up parameter values are set by using test data from the chemical plant or are manually specified.

[00106] As an example, the equations are expressed as follows for the case where no independent hold-up test data is available. The inflow mass rate of any gas or liquid species (i) at the sparger node dm in(i)/dt is determined by the pressure difference between the pressure in the gas feed line (solved outside of the block boundary) and the pressure exerted by the hydrostatic head and column overpressure at the inflow node, at some specific point in time t. The pressure and temperature (and any other relevant) states at each outflow block boundary node are obtained by simultaneously solving the mass and energy conservation equations of the adjacent block. The mass exchange of any species i between the bubble-phase and condensed phase partitions at time t in the column block or any sub-block is represented by dm rx(i)/dt. The outflow rate of the gas molecule from the partitioned bubble-phase at the other side of its boundary is then calculated by minimizing the difference between the target hold-up volume fraction (fVg_target) and the calculated gas hold-up volume fraction (fVg), at time t. Expressed in differential form, the mass conservation equation then states that: dm(i)/dt = dm in(i)/dt + dm rx(i)/dt -km(i) (E3) [00107] Which is expressed in words as the rate of change of the mass of species i is equal to the inflow rate of mass species i plus the rate of change of mass exchange of any species i between the bubble-phase and condensed phase partitions at time t, minus a constant times the mass of species i.

[00108] The symbol k denotes the outflow time constant k (in units of 1/s) is determined by minimizing some objective function (F) of the gas volume holdup difference: k = min F(fVg target -fVg) (E4) [00109] Which is expressed in words as constant k is equal to a minimum of an objective function of the gas volume holdup difference between a target gas hold up and the calculated gas hold up state at a specific point in time.

[00110] The integrated gas species masses in the bubble phase are then converted into the total moles (n), using their respective molecular weights M(i): n = Em(i)/M(i) (E5) [00111] and where the total absolute pressure (P) exerted on the bubble phase in the column block or sub-block volume is then obtained from the average pressures at its inflow and outflow block/sub-block boundary nodes. Since the temperature gradients between the bulk partitioned bubble and slurry phases would be expected to be small and to remain relatively constant over time (the bubble-phase residence time is short compared to the condensed phase), the default approach in this example is to not partition the dynamic energy balance between these phases. This means that the bubble phase temperature (T) can also be averaged for the specific block or sub-block, unless there is strong experimental evidence to support a higher fidelity approach. However, averaging of the P and T state properties in a sub-block becomes less significant as its physical size is reduced. Since the column block is also assumed to be well mixed (maximum mixedness), the simplest case is used as an example of where these average block states are used to calculate the total bubble-phase volume (Vg, at a specific point in time) from the ideal gas law (IGL), with Rg being the ideal gas constant: Vg = n/P *Rg * T (E6) [00112] which is expressed in words as the total bubble-phase volume is equal to the total moles n divided by the pressure, multiplied by the ideal gas constant and the temperature.

[00113] and the average bubble-phase density follows: Sg = Em(i)Ng (E7) [00114] which is expressed in words as the average bubble-phase density equals the sum of the mass of species i divided by the average bubble hold-up volume.

[0011 5] The average gas hold-up in the block or sub-block is then: fVg = Vg/(Vs1 + Vg) (E8) [00116] Which is expressed in words as the average gas hold-up (typically expressed as a fraction or percentage) is equal to the total bubble phase volume divided by the sum of the slurry phase volume and the average bubble density. The slurry phase volume (Vsl) is obtained from the block slurry partition states.

[00117] All of the discussions in this example, up to this point, have been based on first-principles concepts, and where necessary, using simplifying assumptions to average the relevant states over each block or sub-block control volume. The empirical part of this example becomes relevant because of the complex ways individual bubbles behave in the turbulent slurry phase mixtures in the column block, which is rooted in the microscopic scale physics and which requires complex algorithms to describe. However, the above first-principles equations have been formulated so the system behaviour can be modelled in a pragmatic way by using the first-principles mass and energy conservation equations to calculate physics or engineering parameters that are used to design and control process unit blocks. Therefore, these parameters become the link that connects the first-principles equations, describing the macroscopic mass and energy flow behaviour in blocks or sub-blocks, with the phenomenological or empirical equations typically used in the process engineering field. Furthermore, the fact that these first-principles equations are expressed and simultaneously solved in the time domain, means that these parameters/variables become 'describers' of the transient system behaviour and opens up a powerful methodology for high-fidelity modelling and simulation of complex system behaviour.

[00118] The above gas hold-up (fVg) variable is one such example, while the superficial gas velocity (vs) is another variable that is used to correlate observed behaviour (at a macroscopic scale) in bubble columns. The superficial gas velocity can be defined at any point in the block or sub-block where the time-based states are solved explicitly and, again, the smaller the sub-block volumes, the more accurate the described behaviour. In an example case, where the bubble column, as a whole, is treated as the first-principles block, the linear velocity vs is expressed in terms of the average states of the assumed well-mixed inventory material. The numerically (computationally) least expensive option is to assume that most of the physical mass transfer occurs just above the inflow node, that is just above the sparger at the bottom of the column, which then yields: vs = (Zdm(i)/dt)/Sg/As = (Edm(1)/dt)/Zm(i)/Vg/As (E9) [00119] Which is expressed in words as the superficial gas velocity equals a sum of the mass inflow rate of, in this case, gas species i divided by average bubble-phase density, divided by As equals the sum of the mass inflow rate of, in this case, gas species i divided by the sum of the mass of, in this case, gas species i divided by the average bubble-phase volume divided by the cross-sectional area-Where As represents the cross-sectional column area of the bubble swarm, and by combining ES & E6: Vg = Z(m(i)/M(i))/(Ps-Ftg-T) (E10) [00120] Which is expressed in words as the average bubble density equals the sum of the ratio of the mass of species i divided by the molecular mass of species i, divided by the Ps times the ideal gas constant Rg times the temperature. Ps refers to the absolute pressure experienced by the gas bubbles in the region around the sparger (at the bottom of the column) -that is, the atmospheric pressure at that specific geographic location (Patm), added to the block freeboard space overpressure (gas-phase gauge pressure, Pg), added to the column hydrostatic pressure (Ph). This hydrostatics refer to the total combined solid-liquid-gas phase block or sub-block inventory density (Sh = mh/Vh) and its corresponding vertical height (Lh) at that specific point in time, t (g refers to the gravitational constant)-Ps = P atm + Pg + Sh*g*Lh (Ell) [00121] Which is expressed in words as the absolute pressure experienced by the gas bubbles in the region around the sparger is equal to the sum of the atmospheric pressure, the gas-phase gauge pressure and the product of the sub-block inventory density, vertical height and gravitational constant.

[00122] This compression of the bubble-phase volume (Vg) at the column bottom around the sparser, where most of the interfacial mass transfer is assumed to occur, has an important secondary effect in that the partial pressures of all the relevant gaseous molecules (pi) increase in this region (IGL, E6): pi = ni/Vg*Rg*T (E12) [00123] Which is expressed in words as the partial pressures of the gaseous molecules around the sparger is equal to the ratio of the moles of species i to the average bubble volume times the ideal gas constant and temperature.

[00124] Therefore, the thermodynamic states (for example, CO2 eq 8z pH20 eq in El 8z E2, respectively) are higher, as calculated by the embedded first principles and/or data science algorithms (generated from expert thermodynamic software).

[00125] This means that the driving force for interfacial mass transfer is highly transient, and clearly depends on a myriad of interacting effects that are rooted in the material (different species) flows across the block boundaries, the species physical and chemical reaction rates and inventory mass hold-ups within the block boundaries, and the corresponding thermodynamic states (kinetic driving forces of physical & chemical mass transfer). Furthermore, no block is isolated from the environment states like P atm, T atm, etc., and changes in these states at any specific point in time will have casual effects on all of the block steps, including the interfacial mass transfer step (substep 1) discussed here Again, as mentioned above, the smaller the sub-block control volumes, the higher the simulation fidelity.

[00126] An important relationship remaining for substep 1 is to select an appropriate empirical correlation to relate these (or other) variables to the mass-transfer rate at fVg_target. Two examples are given, the first (simplest) case i) representing a low-fidelity, single lumped parameter approach, and the second case, GO representing slightly higher fidelity, which may or may not be justified, depending on what additional information is available for the specific implementation case.

[00127] i) Single lumped parameter example [00128] The physics and chemistry of the fluid surrounding a gas bubble can influence interfacial mass transfer in complicated ways. Various factors are at play, such as the relative fluid/bubble slip velocity, presence of solids particles, liquid ionic strength. The simplest treatment is to use a single lumped parameter that is based on various experimental observations in non-coalescent coarse-bubble heterogeneous flow regime systems, where the gas hold-up and combined mass transfer coefficient both respond as power-law functions to the superficial gas velocity: fVg target = k I -vs0.7 (E13) kLa_target = k2.vs0.7 (E14) [00129] Which is expressed in words as the target gas hold-up is equal to the product of an adjustable parameter kl determined empirically and the superficial velocity, to the power 0.7; and the target average mass-transfer coefficient kLa is equal to the product of an adjustable parameter k2 determined empirically and the superficial velocity, to the power 0.7.

[00130] Since the power value of 0.7 is a typical (relatively insensitive) response observed in non-coalescent coarse-bubble heterogeneous flow systems, kl and k2 are treated as the only adjustable parameters to mimic the results from the specific application in a narrow (known) operating regime where basic mass-transfer type testwork has been conducted to estimate these parameters.

[00131] ii) Higher fidelity example [00132] The simple case above would be limited to a column block or sub-block with very narrow physical and chemical property variations. This second example is an empirical expansion of the above case, with the objective of applying the model to a wider operating range. In its simplest form, the kLa could still be retained as a lumped parameter but then instead of assuming a constant power value of 0.7, we cater for a variable power value (x), depending on the flow regime -that is, varying from the linear functions (x = 1) for the homogeneous flow regime, to varying power laws (x < I, typically 0.85 down to 0.7) for the heterogeneous flow regime, depending on the physics and physico-chemical attributes in the specific control volume (sub-block): n/g_target = kl.vsx (E15) kLa = k2.vsx (E16) [00133] Which is expressed in words as the target gas hold-up is equal to the product of an adjustable parameter k I determined empirically and the superficial velocity, to the power x; and the target average mass-transfer coefficient kLa is equal to the product of an adjustable parameter k2 determined empirically and the superficial velocity, to the power x.

[00134] Where more data is available for a specific apparatus and its contents, higher fidelity can be added. Based on such information, the lumped kLa parameter is then split up, for example, according to the sparger type and coalescent tendency of the bubbles in the specific liquid. If the sub-block control volume is small enough so that the physico-chemical properties remain relatively constant over the specific time step used to solve the first-principles equations (above), but large enough to assume that the equilibrium bubble size is rapidly reached after the gas enters the control volume (or, leaving the sparger,), then deconvolute the kLa into its contributing parts. The total bubble area (a) per unit slurry volume (1/m) should therefore be relatively independent of its spatial position, especially in turbulent mixtures, and can then be related to the average hold-up and mean bubble size (ds) by the following geometrical relation: a = 6/ds. fVg_target /(1-fVg_target) (E17) [00135] Which is expressed in words as the total bubble area is equal to six divided by mean bubble size times the ratio of the target gas hold up divided by one minus the target gas hold up.

[00136] If the experimental hold-up data is supported by average bubble-size measurement, the total bubble area can be obtained, or, an average can be assumed from similar rigid sparger columns and solution ionic strengths, e.g., 5 mm, which then gives an estimate of the total surface area in the control volume block. Similarly, the value of the linear coefficient k.L (m/s) can be calculated from mass transfer experimentation at similar superficial gas velocity and ionic strength solutions, or, can be estimated from the physico-chemical attributes in the control volume, or, simply fixed for coarser average bubble size (>2 mm) in bubble columns with rigid spargers under high mixing intensity, e.g., 4e-4 m/s. In general, the smaller the bubbles, the lower the value of kL. However, the increasing area term (a) typically wins out over the above bubble-size ranges.

[00137] One important mathematical 'connection' remaining for substep 1 in this example, is to connect the interfacial mass transfer to the energy balance of the bulk slurry phase. This is important because the temperature of the interfacial region (Ts), that is, between the bubble phase and the bulk slurry phase, could influence the physical mass transfer rate in non-intuitive ways. Note also that the interfacial temperature (Ts) could be significantly different than the bulk (T) bubble or slurry phase temperatures (or average thereof) (e.g., in E6 above). Therefore, although pi (in E12) could be relatively insensitive to Ts due to the short bubble-phase residence time in the block, the physico-chemical properties in the interfacial layer could be highly sensitive to Ts. For example, decreasing the solubility ([CO2 eq] term in El), increasing the diffusion rate of CO2(a) through the interfacial film layer (i.e., increasing kL), influencing the bubble surface area per unit volume (a) at a microscopic level (ripples) due to changes in surface tension, or even decreasing the total interfacial surface area since bubbles tend to coalesce more at higher temperatures. The mixing intensity of the gas-slurry mixture, the superficial gas velocity and average bubble size would therefore be key aspects influencing the conduction of heat away from the interface. These are clearly complex and highly integrated effects that demand (from a pragmatic perspective) phenomenological and empirical treatment within the mathematical constraints of the first-principles framework.

[00138] Using the lumped parameter example (i) above, the kLa value can then be expressed phenomenologically using the Arrhenius relationship. This is because diffusion processes are also activated, as can be argued from first-principles physical chemistry (e.g., transition-state theory,): kLa(T) = kLa Tr-exp(-Ea/Rg*(14-1/Tr)) (El 8) [00139] Which is expressed in words as the average interfacial mass-transfer coefficient at the block or subblock temperature equals the product of the average interfacial mass-transfer coefficient at the reference temperature and the exponent of negative of the activation energy of molecular diffusion divided by the Ideal Gas Law (gas) constant times the difference between the inverse of the block or subblock temperature and the inverse of the reference temperature.

[00140] Where kLa Tr at the reference temperature (Tr) is known from the lumped parameter approach (E16 -Example (i) above) or calculated using the higher fidelity approach (E17 & related diffusion constants -Example 00 above) at some known bulk-phase temperature (T) in the column block or sub-block. Typical measured activation energies (Ea) for an interfacial gas-liquid diffusion step (substep 1), in the presence of a chemical reaction step (substep 2 -see next section; at the other side of the interface), vary roughly between 15 and 35 kJ/mol.

[00141] Alternatively, an empirical equation, can be used, based on the observation that the Itha typically increases the temperature by about 2%/°C: kLa(T) = kLa Tr-0(T-Tr) (E19) [00142] Which is expressed in words as the average interfacial mass-transfer coefficient at the block or subblock temperature equals the product of the average interfacial mass-transfer coefficient at the reference temperature and the temperature correction factor and the difference between the block or subblock temperature and the reference temperature.

[00143] An example value of the temperature correction factor 0 is 1.024, but values as low as 1.01 and as high as 1.09 have may be used. These two equations (E18 & E19) yield very similar functionalities in the temperature response and can be used to relate the kLa at the bubble surface to the bulk temperature as follows: kLa(Ts) = kLa(T)*0(Ts-T) (E20) [00144] This is an example to show how effects in the environment, such as the diurnal cycle, can be connected to microscopic level phenomena that may have important bearing on the interfacial mass-transfer rate (substep 1) via the block or sub-block bulk temperature (T).

[00145] The above examples give some understanding of the different methodologies that are followed to add fidelity to the model for the gas-liquid mass-transfer step. Importantly, it shows how the complex physical and chemical phenomena can be embedded in the first-principles block framework using phenomenological and empirical correlations -that is, to capture the complex physics and chemistry within the algebraic constraints imposed by the transient first-principles mass and energy flows. This unique methodology also allows the calculation of physical and chemical properties, and related lumped parameters (typically used by these correlations), at points in the column block where it would be difficult to measure experimentally. Notwithstanding, the more experimental data available (wherever measurement is indeed practically possible), the higher the modelling fidelity and the more algebraic constraints imposed by this methodology, and hence, the more predictive the behaviour of the interfacial mass-transfer step (substep 1). The next two steps (below), namely the chemical aqueous sodium bicarbonate NaHCO3(a) formation (substep 2) and solid sodium bicarbonate NaHCO3(s) precipitation (substep 3), follow the same methodology. Examples are however used to demonstrate that all three steps are coupled physically through various acausal phenomena which are accounted for mathematically.

[00146] ,S7oichiornetric chemical reaction of salts (aqueous reaction substep-2) [00147] Substep 2 connects physically to substep 1 via the CO2(a) concentration at the liquid-side of the bubble phase interface, and any other slurry surface gas/liquid interfaces, and is represented by the following overall chemical reaction stoichiometry: CO2(a) + H20(1) + Na2CO3(a) <=> 2NaHCO3(a) (R3) [00148] Which is expressed in words as aqueous carbon dioxide plus liquid phase water plus aqueous sodium carbonate is in equilibrium with aqueous sodium bicarbonate.

[00149] The physical transfer rate of CO2(a) across the gas-liquid interface could therefore be an important resistance (limiting factor) to this chemical (R3) reaction rate. But, the relative importance of this resistance depends on the chemical potential (at a quantum level) of R3 to move in a forward (left to right) direction. When looking at this from the reverse (right to left) direction, the precipitation rate of NaHCO3(a) => NaHCO3(s) (substep 3) will also impart a resistance on the net forward rate, because the slower this precipitation rate, the higher the aqueous NaHCO3(a) concentration and the lower the chemical potential of R3 to move in a net forward direction at an economical rate, at any point in time (t).

[00150] To mathematically embed this 'reversible' chemical step (substep 2) into the time-based (dynamic) first principles mass, energy and particle population conservation framework, requires understanding of what actually happens at a molecular level. A more rigid understanding of the molecular interactions is developed by separating the above stoichiometric reaction (R3) into its forward and backward directions: CO2(a) + H20(1) + Na2CO3(a) => 2NaHCO3(a) (R3a) 2NaHCO3(a) => CO2(a) + H20(1) + Na2CO3(a) (R3b) [00151] The overall stoichiometry is still accounted for at the macroscopic (measured) scale, but the varying state properties of the bulk liquid phase, as described in this transient column block or sub-block example, will cause these individual (R3a, R3b) reaction rates to change in complex ways. This means that by using the overall equilibrium reaction (R3) representation as the departure point for mathematically describing its net forward rate, the acausal physical-chemical interactions at the liquid-side of the bubble phase (substep 1) and the acausal chemical-physical interactions at the liquid-side of the particle phase (substep 3) are not taken into account. It is therefore clear that only if these chemical interactions are i) slow enough to occur in the bulk control volume, that is, physically (far enough) removed from the gas-liquid (Step I) and liquid-solid (Step 3) interfaces, and, ii) fast enough so that the net reaction R3 is at pseudo-equilibrium (i.e., Rate of R3a Rate of R3b), then the bulk species concentrations (over some small numeric time step At, at time t) can be used to describe the net forward chemical reaction rate of R3 via transition-state theory along the reaction coordinate that is valid for that specific time step. In its simplest form, the R3 rate can then be treated phenomenologically as a reversible homogeneous reaction, but without making any molecularity assumptions, that is -the reaction orders (n) are assumed to be unrelated to the reaction stoichiometry and hence viewed as non-elementary reactions. Furthermore, if it is assumed that the temperature dependencies of the pre-exponential constants (k) and the activation energies (Ea) are small compared to the temperature term in the exponent, the same functional forms are obtained, irrespective of whether the first-principles basis is transition-state or collision theory. This then leads to the following rate expression for R3: k(i) * dm rx(i)/dt = kf* exp(-Ea f/Rg*(1/T-1/Tr)) -kb * exp(-Ea b/Rg*(1/T-1/Tr)) (E21) [00152] where species (i) refers to any of the formal bulk-phase species CO2(a), H20(1), Na2CO3(a) or NaHCO3(a) taking part in the reaction, and where k(i) is the proportionality constant for the specific species, so that the rate of its bulk-phase mass change (due to the chemical reaction) in the column block or sub-block can be related to the overall reaction (R3) stoichiometry (molar basis). The rate constants for the forward and backward directions over time step At (at time t) at the reference temperature (Tr) are then, respectively: kf = kl * [CO2(a)]111 * [Na2CO3(a)]n2 (E22) kb = k2. [NaHCO3(a)]n3 (E23) [00153] where the square brackets again refer to the time-integrated species concentrations in the bulk liquid. If enough data is available over carefully controlled conditions in the laboratory, the 7 unknowns (Ea f, kl, nl, n2, Ea b, k2, n3) can be obtained via regression However, even then, after careful parameterization of this phenomenological rate equation (E21), its applicability to the multi-phase transient control block mass and energy conservation framework is questionable -the reason being, the many assumptions required to derive E21, E22 and E23, in addition to their narrow application limits when neglecting the acausal physical-chemical interactions.

[00154] At the other modelling extreme, the premise can be used that the (chemical) potential energy is the most fundamental (quantum-level) driving force description of solvated species interactions, and which is consistent with a thermodynamic treatment, using ab initio modelling and the reference state of an electron at rest in vacuo). That is, it assumes that the kinetic energy (translational, rotational, etc.) cancels approximately between the solvated reactant and product species at constant temperature. It can therefore be argued that the driving force (rate) for R3 (whether moving from left to right, R3a, or vice versa, R3b) can be expressed as the Gibbs-energy based deviation from the energetically ideal equilibrium state. Since the time to reach this thermodynamic reference state is significantly shorter than the mass and energy changes of the bulk (block or sub-block) control volume, it does not warrant dynamic (time-derivative) treatment within this block mass, energy and particle population conservation framework, besides the fact that it would be computationally cumbersome (numerically stiff, due to the significant difference in the time scales of the macroscopic flow algorithms and these microscopic rate phenomena). This thermodynamic reference approach, to describe the chemical rate (Step 2), is akin to using the thermodynamic solubility state to describe the rate of interfacial (physical) gas molecule mass transfer (substep I) and using the thermodynamic supersaturation state to describe the rate of precipitation (substep 3). As mentioned for substep 1 (previous section), these thermodynamic state calculations are ideally suited for treatment using specialized solution-thermodynamic software and machine learning tools to generate and capture that information separately, and then to embed it as a stand-alone data-science block into the first-principles framework via pre-defined input/output (I/O) ports.

[00155] The phenomenological reaction-kinetic approach (E21, E22 and E23) can be combined using the potential energy argument to reduce the number of adjustable parameters. Although this approach is largely empirical, it allows early evaluation of pilot/demo plant results over narrow operating regimes with less computational effort. The rate constant is then expressed as a combined potential energy term via some function ft[CO2(a)], [Na2CO3(a)], [NaHCO3(a)]), and using the assumption that no significant concentration gradients develop between the bulk liquid and the gas-liquid interface. The influence of the 'reverse' rate is then implicitly taken into account by expressing the empirical rate constant kc (1/s), at time t, as: kc = kr.([1\la2CO3]/NaHCO3(a)Dn.exp(-Ea/Rg*(1/T-1/Tr)) (E24) [00156] and the phenomenological rate expression then becomes: dn rx/dt = kc.[CO2(a)] (E25) [00157] This first-order dependence on the dissolved CO2 concentration is used to link the bulk chemical reaction rate to the dynamic gas-liquid interface phenomena in the bubble phase. The substep 1 mass transfer equation is now briefly revisited, but without assuming that the above chemical reaction (substep 2) is spatially separate from the diffusion layer (substep I), which means that their respective algorithms cannot be solved independently.

[00158] The reality (demonstrated both theoretically and experimentally) is that when the chemical reaction (substep 2) is fast relative to physical diffusion (substep I), the chemical potential energy (quantum-level driving force) accelerates the diffusion rate because the chemical reaction front moves inwards (protrudes into the diffusion film, making it thinner'). Also, when the chemical rate does not dominate the diffusion rate (e.g., in the absence of a catalyst), then aqueous species concentration gradients would be expected to develop at the liquid-side of the interface. Although merged (physical/chemical) rate expressions can be developed for theoretically ideal cases the real-world physics at a microscopic (molecular chemical reaction & diffusion) level and at a macroscopic (transient liquid and gas flows) level create complex scenarios that are not static. However, the powerful mass, energy and particle conservation algorithms used to describe the physical behavior in this example, lays down the mathematical platform needed to digitally describe these complex and transient phenomena at high fidelity.

[00159] Again, the data-science approach provides the most accurate physico-chemical attributes and chemical potential expression via the calculated Gibbs free energy, either using the bulk mass and energy states, or, using the calculated or inferred states at the gas-liquid interface. Expressed differently, these first-principles frameworks provide the algebraic constraints needed so that each data-science block can be trained using narrower, well-defined input parameter spaces. Either way, whether using the advanced data-science approach, or, instead, an initial-phase (low-fidelity) empirical approach, such as E24, the chemical reaction rate equation (such as E25) is mathematically connected with the diffusion rate equations (El 8-E20, etc.). The higher the fidelity of the molecular chemistry description (e.g., by using the data-science implementation), the simpler (less parametrized) the phenomenological/empirical algorithms required to link the gas-liquid interface behavior with the macroscopic (operational) parameters, such as the superficial gas velocity (vs), and with the microscopic phenomena, such as the diffusion-chemical reaction dependencies described above.

[00160] Using the simplest phenomenological approach (E24-E25), as an example case, the link between the chemical and diffusion rate expressions is made by defining an enhancement factor Ec by which the chemical reaction increases the rate of CO2(g) absorption, as compared to absorption without any reaction. With evidence that the R3 chemical reaction responds proportionally (first-order) to the change in [CO2(a)], the dimensionless enhancement factor is then: Ec = (I + kc.D/kL2)0.5 (E26) [00161] where D is the CO2(a) molecular diffusion coefficient (m2/s). The CO2(g) absorption rate then becomes: dn(CO2)/dt = Ec*kLa(T).([CO2 eq] -[CO2(a)]) (E27) [00162] In their simplest forms, these equations refer to the bulk temperature T and aqueous species concentrations. This means, that the rates of thermal and species diffusivities between the partitioned bubble and the bulk phases in the block or sub-block are fast. However, such gradients (e.g., Ts), when supported by physical (real-world) observations, are easily incorporated using the above framework, and based on the net rate of local heat (enthalpy) addition (due to the absorption, E27) minus the thermal diffusivity rate to the bulk phase.

[00163] Importantly, the above equations example provides the means to connect the dynamic changes in the block or sub-block environment with its transient chemical behavior (substep 2), which, in turn, is coupled to the acausal mass and energy flow behavior of the partitioned bubble and slurry phases, as well as the transient chemico-physical interactions at the gas-liquid interface (substeps I & 2).

[00164] Precipitation from liquid to solid (precipitation of salts substep-3,) [00165] With the required mathematical descriptions in place for the block interfacial mass and energy transfer rates (substep I), and, the associated homogeneous chemical reaction mass and energy rates (substep 2), the only part remaining is to account for the solids-phase behavior. That is, so that the first-principles column block or sub-block mass and energy conservation framework can dynamically connect the macroscopic mass and energy states to the microscopic (quantum-level) driving forces and molecular-level mass and energy changes. The first-principles network of algorithms, describing the column block or sub-block behavior in terms of measurable state properties, is then complete so that the more complex microscopic phenomena can be treated empirically or phenomenologically, or where warranted, using machine learning in the data-science domain.

[00166] Substep 3 connects (mathematically) to substep 2 via the NaHCO3(a) concentration at the liquid-side of the liquid-solids interface. It is represented by the following overall chemical reaction stoichiometry: NaHCO3(a) <=>NaHCO3(s) (R4) [00167] Again, as was the case with the gas phase behavior at the gas-liquid interface (Steps 1-2), the rate (kinetics) of NaHCO3(s) surface formation (Steps 2-3) is built up from first-principles interfacial energy concepts. That is, as some function of the relative deviation from thermodynamic equilibrium. But before expressing the microscopic driving forces behind salt precipitation and dissolution, the macroscopic-level numerics are added to complete the first-principles algorithm network in terms of the bulk mass/energy transient flows across and within the block. But unlike the previous discussions, the substep 3 precipitation reaction rate requires transient (time domain) treatment in the particle volume dimension. The reason for this lies in the energetics. More work is required to form new particles (nucleation) as compared with forming larger particles (growth). These transient particle-level rate phenomena influence the mass and energy flows at the liquid-solid interface (substep 2-3), and therefore, also across the whole block or sub-block first-principles algorithm network.

[00168] The theoretical background is based on the mass conservation along the characteristic particle size distribution (PSD) coordinate x. The particle size density distribution Dx is then: dN = Dx dx (E28) [00169] where N represents the total number of particles in the block or sub-block control volume The one-dimensional mass conservation (continuity) equation is used to close the particle balance: axiat + a(Dx-v)/ax = aDx.in/at -aDx.out/at + aDx.13(x)/at -aDx.1)(x)/6t (E29) [00170] The in-and outflow terms, Dx,in and Dx,out, respectively, account for the PSD changes due to the slurry entering and exiting the block or sub-block bulk control volume. Since the gas absorption (substep 1), the chemical reaction (substep 2) and the precipitation (substep 3) generate heat, and therefore temporarily raise the local temperature (before dissipating into the bulk material via convection and conduction), it can be assumed that the substep 1-2 coupling would dominate in the partitioned bubble phase, whereas the substep 2-3 coupling would tend to dominate in the bulk slurry phase, especially in the cooler areas closer to the column wall. The above mass balance (E29) is therefore defined per unit mass bulk slurry, which excludes the bubble phase.

[00171] The internal velocity term (v) describes the rate (m/s) of particle growth (or shrinkage) due to crystallization (or re-dissolution) in the column block or sub-block slurry phase. The net particle birth (Dx,B) and death (Dx,D) contributions to the density distribution only vary with time (t) because their respective rates are dominated by the (integrated) particle mass in the discretized size coordinate. The reason for this is rooted in the probability that particles of a specific size class (Ax) will meet and stick together with particles of the same or other size classes. The same can be said for particles of a specific size class (Ax) bumping into other particles and objects to form smaller particles. That is, particulate mass breaks from larger to smaller size classes. These birth and death rate contributions to the transient particle density changes cannot yet be formulated from first principles theory because these phenomena are brought about by localized fluid dynamics and kinetics. However, the way in which the column block transient mass and energy changes are set up, using the above first-principles network of algorithms, allows sub-block partitioning for the sole purpose of quantifying those phenomena in a more predictable manner.

[00172] An example of this would be where the sub-block control volume in the column is defined to represent the gas-free bulk slurry of specified thickness in contact with the column wail. This sub-block would then experience slightly lower temperature at night, i.e., higher NaHCO3(s) supersaturation, due to its energy (convective and conductive heat) transfer to the environment. In addition, the less turbulent and lower fluid velocity at the wall would then increase the probability of agglomerate formation. For this example, the assumption is made that the agglomerates represent collections of weakly bound particles where the resulting mass is similar to the sum of the individual crystal masses. This means that insignificant growth is associated with the formation of agglomerates comprising weakly fused crystals.

[00173] The birth rate term, 5Dx,B(x)/Dt, in E29 can therefore be split into (i) a source rate term, due to nucleation at the minimum particle size x(1), 0Dx,B(x(1))/Ot = J, and, (ii) an aggregation kernel, kercaDx,agg(x)/00 = A(x), i.e., the population density distribution changes with x such that the net mass change due to agglomeration remains zero over time t. An analogous breakage kernel is used to express the death rate term (in E29) due to crystal breakage and attrition, ker(aDx,B(x)/a) = B(x), which, again, preserves the mass balance over the x coordinate as particles break (probabilistically) out of larger size classes into smaller size classes. Both the aggregation and breakage can therefore be expressed explicitly in the size (x) coordinate: (Mx/at + a(Dx-v)/ax = anx,in/at -aDx.outiat + J + A(x) + B(x) (E30) [00174] where J essentially becomes a boundary condition at the smallest particle size class when numerically solving this equation in the discretized x coordinate. The velocity term (v), on the other hand, cannot be assumed to be constant over size. The geometry/size complexities arise due to the way polycrystalline aggregate structures, e.g., the remnants of the agglomeration, are often accompanied by an increase in particle mass. That is, because if the attractive forces prevail long enough, crystalline bridges (due to growth) will cement the particles together. This phenomenon could become more prevalent at high supersaturations so that the inter-particle cemented aggregates become more compact and remain intact, even under high slurry turbulence. A phenomenological description can be formulated by lumping these geometric/mechanistic complexities together as a one-dimensional effect, which means that the velocity term is re-casted into the one-dimensional characteristic size along the x coordinate: v = vx = dx/dt = G*xn (E31) [00175] This is essentially a topochemical reaction boundary assumption lS2012al, and, although it imposes an idealized surface area on the particle size domain, it helps to account for these complexities when incorporated as some (empirical) function of the supersaturation S, i.e., n = f(S), which could, for example, vary between 0 (zero order) and 1 (first order) in the specific column block or sub-block slurry phase control volume.

[00176] The above formulations (E30 & E31) divide the growth of NaHCO3(s) (via R4) into two distinct theoretical mechanisms: (i) heterogeneous nucleation and rapid associated growth (J), which then dissipates (as rapidly) when the particle size of the smallest class Ax(1) is reached, and, (ii) crystal growth on existing NaHCO3(s) particles surfaces, which is numerically dealt via particle fluxes that automatically conserve the number of particles (population balance). There are a myriad of nucleation and growth mechanisms, and their specific functional forms, that can be implemented into the above conservation equation. However, the transient nature of the fluid dynamics, and varying temperature, chemical species and particle surface area in the column block or specific sub-block would likely not accommodate such highly specific mechanistic descriptions, unless the operating state of the block can be carefully controlled. Hence the above default approach that assumes 'instantaneous' growth associated with the nucleation number density term (J) in the smallest particle size class, where 'instantaneous' means such (limited) growth is much faster compared to the time scale of the crystal growth velocity term (G). Also, if the smallest size class, Ax(1), is defined small enough, then the solids mass change associated with implicitly lumping the rapid growth to the J term becomes insignificant. The main purpose of the nucleation term J is therefore to account for the addition of new particles to the number density distribution, whereas the growth rate term G (m/s) focuses on the way the existing particles grow larger and account for the associated mass change (slurry density) of the solids phase.

[00177] These two micro-model terms are however closely related in the way they affect the PSD and crystallization rate in the column block or sub-block. The supersaturation, S. which describes how 'far' the system is from thermodynamic equilibrium, is the most important variable that determines the relative rates of the nucleation and growth terms. Therefore, as with the gas-liquid mass transfer (Step 1) and chemical (Step 2) rates, the thermodynamic driving force (here expressed as 5) determines the kinetics (rate) of interfacial liquid-solid mass transfer: S = [NaHCO3(a)]/[1\1aHCO3 eq] (E32) [00178] Again, as with substep 1 and 2, these thermodynamic state calculations are ideally suited for treatment outside the above framework -that is, using specialized solution-thermodynamic software and machine learning tools to generate and capture that information separately, and then to embed it as a stand-alone data-science block into the first-principles framework via pre-defined input/output (I/O) ports.

[00179] The dependencies of the nucleation J and growth G rate terms on the supersaturation are related to the mechanism on a molecular level. This follows directly from basic crystallization theory. The J and G rates typically exhibit high and lower-order dependencies on S, respectively. Since there are many different nucleation and growth mechanisms, with complex dependencies in S, simple phenomenological-type descriptions are most appropriate for this example. Firstly, both particle nucleation and crystal growth algorithms are fundamentally rooted in the free energy changes as the solute moves between the liquid and solid state. These are activated processes related to the probability of crossing the activation energy barrier and, again, like with substep 1 and 2, the well-known Arrhenius equation is used to account for the temperature effect. Since it takes more energy to create new solids surface via nucleation, cf growing on existing surfaces, the order (n) of the dependence on S is typically higher for.1-and lower for G. The following two equations are semi-empirical examples to describe these terms: J = kj * Sn. exp(-Ea/Rg*( 1/T-1/Tr)) (E33) [00180] where kj is the rate constant (1/s), expressed as a number density contribution to the smallest size class in the PSD, with a typical value for n = 4.

G = kx * ( S-1)n * exp(-Ea/Rg*(1/T-1/Tr)) (E34) [00181] where kx is the linear growth rate constant (m/s) and S-1 is the relative supersaturation, with a typical value for n = 1 or 2.

[00182] The final two phenomena to describe are the aggregation A(x) and breakage B(x) terms. As pointed out above, these are probabilistic contributions to the transient particle density changes, which are difficult to formulate from first principles theory. However, since their usage is defined within the first-principles network of block or sub-block algorithms and the number balance framework described above, these terms can be treated empirically or quantified using data science tools, such as machine learning. The general descriptions below are examples of how these terms could be parameterized: A(x) = ka*Pa(x)*Da(x) (E35) [00183] where ka is the rate constant (Us), Pa is a probability vector (at time t over the x coordinate), based on some weighting factor times the fraction pulp density of each size class at time t. The distribution vector, Da, describes the change in the number density distribution (over the x coordinate) so that the total particle mass (volume) is conserved -that is, as the smaller particles agglomerate into larger sizes according to some ideal, e.g., spherical, geometry. The analogous equation for the death rate (breakage/attrition) term is then: B(x) = kb*Pb(x)*Db(x) (E36) [00184] where kb is the rate constant (1/s), Pb is a probability vector, and Db is the number density distribution vector, with opposite sign (direction) to Da because larger particles break into smaller fragments. The weighting factors of the probability vectors (Pa and Pb) are best suited for integration with machine learning tools, that can be trained using PSD measurements in specific block or sub-block control volumes, e.g., where diurnal cycles in the environment temperature tend to increase agglomeration, or, where the turbulence or other forces in the specific block or sub-block are large enough to break larger aggregates or particles into smaller fragments. But, importantly, these machine learning blocks are trained and implemented within the above first-principles network of algorithms, which helps to account for these types of microscopic complexities in a simpler and more predictable manner. Also, looking at this from a macroscopic process perspective, the better the overall integration of the different blocks are controlled, the more accurate and predictive these data-science blocks become.

[00185] The above algorithm framework links the substep 3 solids formation rate with the substep 2 chemical reaction rate, based on the thermodynamic driving force (via supersaturation S, at temperature T), the particle number conservation (via solving for 0Dx/at) and the block or sub-block states (via its bulk mass and energy balance conservation equations), at each time step, t. The complex substep 2-3 rate phenomena are therefore just as complex as the substep 1-2 rate phenomena. This means that the substep 1-2-3 phenomena are highly integrated, often making the system respond counter-intuitively, including the way it interacts with the environment during the diurnal cycle. A tangible impact on the efficient SBC column control can therefore made by using the above algorithm network and data (literature, databases, laboratory, pilot, demo and full scale) to understand and improve each site-specific design, implementation and control, related R&D, and techno-economic evaluation over a complex and varying parameter space, including quantifying and predicting how the site-specific environment, such as the temperature and wind speed, would affect each operation.

[00186] In FIGs. 11 to 22 the following notation is used: First principles FPR Phenomenological PHN Empirical (entirely) EMP Machine learning MLN Numeric error minimization NEM Time T Total pressure P Temperature T Total mass M Species (i) Mass m Species (i) Partial Pressure p Total enthalpy H Species (i) Enthalpy h Particle number density distribution D Difference A Total derivative d Partial derivative d Block states following integration Dark filled rectangle with rounded corners PHN, EMP, NEM algorithms, FPR algorithms Medium filled rectangle with rounded corners Algorithm, PRIM, EMP operation or assignment Unfilled rectangle with rounded corners Comment or description Unfilled rectangle without rounded corners Dashed arrow Causal or acausal flow (visual aid) Solid arrow Rate (variable) assignment Dot dash arrow State variable assignment Dotted arrow Comment or variable (visual aid) Solid small circle Equation or variable reference point (indicative visual aid) Large white circle Gas bubble (indicative visual aid) Annulus Gas bubble with diff film layer (indicative visual aid) Solid medium circle Solid particle (indicative visual aid) Double headed wide arrow Chemical or physical reaction (indicative visual aid) Diamond Signal splitting/confluence Filled small rectangle In-/outflow node Unfilled small rectangle Input signal port E Equation mentioned in the specification, where the number following the letter E denotes the number of the equation Node A Mass & energy conserving port of the column example where the gas phase enters the block Node B Mass & energy conserving port of the column example where the slurry phase enters the block Node C Mass & energy conserving port of the column example where the gas phase exits the block Node D Mass & energy conserving port of the column example where the slurry phase exits the block [00187] FIG. 11 is an example of a first principles framework of Block-1 of the digital twin of FIG. 9. Regions of the framework are numbered as follows 1000, 1002, 1004, 1006, 1008. Section 1000 is shown enlarged in FIG. 12, Section 1002 is shown enlarged in FIG. 13. Section 1004 is shown enlarged in FIG. 14. Section 1006 is shown enlarged in FIG. 15. Section 1008 is shown enlarged in FIG. 16.

[00188] FIGs 12 to 16 are enlarged portions of FIG. 11 as explained above.

[00189] The reactor vessel is represented as a rectangle 1010 spanning figures 1000 and 1004. Within the reactor vessel bubbles rise up from the base of the reactor vessel and this is represented using large unfilled circles in FIGs. 11, 12 and 14. In practice there are many bubbles although only four are shown for clarity. A film layer around each bubble is modelled as explained above and the film layer is depicted in FIGs. 11, 12 and 14 as an annulus around the large unfilled circles. Liquid slurry is in the reactor vessel and is represented below line 1012. Above line 1012 is gas.

[00190] From FIGs 11 to 16 it is seen that there are three mass balance equations, one particle population balance equation and one energy balance equation as listed below. These are referred to as conservation equations. Solutions to the phenomenological, empirical and machine learning relations of the corresponding framework in FIG. 17 are constrained by the conservation equations. The use of the conservation equations indicated in FIGs. 11 to 16 is found to be particularly effective.

[00191] Gas sub block mass balance conservation equation is: dm/dt = Zdm/dt(in) -Zdm/dt(out) + Zdm/dt(rxn) [00192] which is expressed in words as the rate of change of mass in the gas sub block over time is equal to the sum over all species i involved in the process, of the rate of inflow of mass of species (i) minus the sum over all species i of the rate of outflow of mass of species (i) plus the sum of the rate of change of mass of species i due to the subblock physical and chemical reactions (rxn) [00193] The bubble phase mass balance conservation equation is: dm/dt = Zdm/dt(in) -Edm/dt(out) + Zdm/dt(rxn) [00194] Thus the rate of change of mass of the bubble phase is expressed in the same way as for the gas sub block mass balance.

[00195] The particle population balance conservation equation is: anxiat + a(Dx-v)/ax = aDx(in)/at -aDx(out)/at + J + A(x) + Mx) [00196] Which is expressed in words as the particle population density distribution (along the internal, one-dimensional, particle size coordinate) partial derivative (change) with respect to time in the block or subblock plus the partial population density derivative of the product of the particle density distribution and the growth rate (along the internal, one-dimensional, particle size coordinate as a complete time derivative), with respect to the particle size coordinate equals the partial derivative of the particle population density distribution with time due to the inflow of solid particles minus the partial derivative of the particle population density distribution with time due to the outflow of solid particles plus the contribution, J, to the particle population density distribution due to formation rate (as a complete time derivative) of new particles of fixed (smallest) size via nucleation plus an aggregation kernel, A, that takes care of the particle population density distribution changes in x (as a complete time derivative) such that the net mass change rate due to agglomeration remains zero over time plus a breakage kernel, B, that takes care of the partial particle population density distribution changes in x (as a complete time derivative) such that the net mass change rate due to breakage remains zero over time.

[00197] The column block energy balance conservation equation is Em-cp(T) dT/dt -dQ/dt = ZdH/dt(in) -ZdH/dt(out) + Edm/dt-h(T) [00198] Which is expressed in words as the product of the sum of the product of each species i mass and their specific heat capacity at the block or subblock temperature and the complete derivative of the block or subblock temperature with respect to time minus the change (derivate with respect to time) of heat exchange of the block or subblock with its environment equals the total enthalpy change due to the sum of all the block or subblock inflow streams minus the total enthalpy change due to the sum of all the block or subblock outflow streams plus the sum of the product of each species i mass change (complete derivate with respect to time) and their specific enthalpy at the block or subblock temperature.

[00199] The slurry phase mass balance conservation equation is dm/dt = Zdm/dt(in) -Edm/dt(out) + Zdm/dt(rxn) [00200] Thus, the rate of change of mass of the slurry phase is expressed in the same way as for the gas sub block mass balance.

[00201] In various examples, at least one of the blocks of the digital twin represents process in a reaction vessel in the chemical plant, and the mass and energy conservation equations for the at least one block are: a gas subblock mass balance conservation equation which is expressed in words as the rate of change of mass in a gas subblock of the block over time is equal to the sum over species i involved in the process, of the rate of inflow of mass of species (i) minus the sum over species i of the rate of outflow of mass of species (i) plus the sum of the rate of change of mass of species i due to subblock physical and chemical reactions; a bubble phase mass balance conservation equation expressed in the same way as for the gas subblock mass balance conservation equation, a particle population balance conservation equation which is expressed in words as a particle population density distribution partial derivative with respect to time in the Nock or subblock plus a partial population density derivative of a product of the particle population density distribution and a growth rate, with respect to a particle size coordinate equals a partial derivative of the particle population density distribution with time due to inflow of solid particles minus a partial derivative of the particle population density distribution with time due to outflow of solid particles plus a contribution to the particle population density distribution due to formation rate of new particles via nucleation plus an aggregation kernel that takes allows for the particle population density distribution changes such that a net mass change rate due to agglomeration remains zero over time plus a breakage kernel that takes care of partial particle population density distribution changes such that a net mass change rate due to breakage remains zero over time; a column block energy balance conservation equation which is expressed in words as a product of the sum of the product of each species i mass and specific heat capacity at the block or subblock temperature and the complete derivative of the block or subblock temperature with respect to time minus the change of heat exchange of the block or subblock with its environment equals the total enthalpy change due to the sum of the block or subblock inflow streams minus the total enthalpy change due to the sum of the block or subblock outflow streams plus the sum of the product of each species i mass change and their specific enthalpy at the block or subblock temperature; and [00202] a slurry phase mass balance conservation equation expressed in the same way as for the gas subblock mass balance conservation equation. Using this particular combination and specification of conservation equations is found to give particularly good fidelity.

[00203] FIG. 17 is an example of a phenomenological, empirical and machine learning framework corresponding to the first principles framework of FIG. 11 and for use in a digital twin such as that of FIG. I. The use of the conservation equations indicated in FIGs. 11 to 16 is found to be particularly accurate and efficient when used with the framework of FIG. 17.

[00204] FIGs 18 to 22 are enlarged portions of FIG. 17. Regions of the framework are numbered 2000, 2002, 2004, 2006, 2008. Section 2000 is shown enlarged in FIG. 18. Section 2002 is shown enlarged in FIG. 19. Section 2004 is shown enlarged in FIG. 20. Section 2006 is shown enlarged in FIG. 21. Section 2008 is shown enlarged in FIG. 22.

[00205] In the example framework of FIGs. 17 to 22 the following equations, which have been explained earlier in this document, are used: El, E2, E4, E6, E8, Ell, E15, E16, E17, E18, E19, E24, E25, E26, E27, E31, E32, E33, E34, E35, E36. Thus, a large number of equations are solved simultaneously for a given time step and are constrained to agree with the conservation equations of FIG. 11 to 17. In this way an efficient and accurate simulation is achieved.

[00206] FIG. 23 shows a computing-based device in which a digital twin is implemented.

[00207] In various examples there is a digital twin for use in controlling an automated chemical plant, the digital twin comprising a plurality of connected blocks, each block representing part of a process carried out by the automated chemical plant; for each block, mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block; for each block, one or more equations describing the part of the process in the block and formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships; and wherein the digital twin is configured such that, for each block, outputs of the equations are constrained to agree with the mass and energy conservation equations; and wherein the digital twin has an output for sending simulation results to a SCADA for controlling the chemical plant using results of the simulation.

[00208] The use of mass and energy conservation equations as well as second equations in the digital twin enables the control of the chemical plant in an unconventional manner to achieve production efficiencies.

[00209] The digital twin improves the functioning of the underlying computing device for control of the chemical plant by use of mass and energy equations in combination with second equations.

[00210] FIG. 23 illustrates various components of an exemplary computing-based device 3008 which are implemented as any form of a computing and/or electronic device, and in which embodiments of a digital twin 104 are implemented in some examples.

[002111 Computing-based device 3008 comprises one or more processors 3000 which are microprocessors, controllers or any other suitable type of processors for processing computer executable instructions to control the operation of the device in order to simulate a process in a chemical plant and optionally to control the chemical plant using the simulation. In some examples, for example where a system on a chip architecture is used, the processors 3000 include one or more fixed function blocks (also referred to as accelerators) which implement a part of the method of any of FIGs. 2 to 22 in hardware (rather than software or firmware) Data store 3004 holds parameter values, equations, measurements and other data.

[00212] The computer executable instructions are provided using any computer-readable media that is accessible by computing based device 3008. Computer-readable media includes, for example, computer storage media such as memory 3002 and communications media. Computer storage media, such as memory 3002, includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or the like. Computer storage media includes, but is not limited to, random access memory (RAM), read only memory (ROM), erasable programmable read only memory (EPROM), electronic erasable programmable read only memory (EEPROM), flash memory or other memory technology, compact disc read only memory (CD-ROM), digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transmission medium that is used to store information for access by a computing device. In contrast, communication media embody computer readable instructions, data structures, program modules, or the like in a modulated data signal, such as a carrier wave, or other transport mechanism. As defined herein, computer storage media does not include communication media. Therefore, a computer storage medium should not be interpreted to be a propagating signal per se. Although the computer storage media (memory 3002) is shown within the computing-based device 3008 it will be appreciated that the storage is, in some examples, distributed or located remotely and accessed via a network or other communication link (e.g. using communication interface 3006).

[00213] The term 'computer' or 'computing-based device' is used herein to refer to any device with processing capability such that it executes instructions. Those skilled in the art will realize that such processing capabilities are incorporated into many different devices and therefore the terms 'computer' and 'computing-based device' each include personal computers (PCs), servers, mobile telephones (including smart phones), tablet computers, set-top boxes, media players, games consoles, personal digital assistants, wearable computers, and many other devices.

[00214] The methods described herein are performed, in some examples, by software in machine readable form on a tangible storage medium e.g. in the form of a computer program comprising computer program code means adapted to perform all the operations of one or more of the methods described herein when the program is run on a computer and where the computer program may be embodied on a computer readable medium. The software is suitable for execution on a parallel processor or a serial processor such that the method operations may be carried out in any suitable order, or simultaneously.

[00215] Those skilled in the art will realize that storage devices utilized to store program instructions are optionally distributed across a network. For example, a remote computer is able to store an example of the process described as software. A local or terminal computer is able to access the remote computer and download a part or all of the software to run the program. Alternatively, the local computer may download pieces of the software as needed or execute some software instructions at the local terminal and some at the remote computer (or computer network). Those skilled in the art will also realize that by utilizing conventional techniques known to those skilled in the art that all, or a portion of the software instructions may be carried out by a dedicated circuit, such as a digital signal processor (DSP), programmable logic array, or the like.

[00216] Any range or device value given herein may be extended or altered without losing the effect sought, as will be apparent to the skilled person.

[00217] Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims.

[00218] It will be understood that the benefits and advantages described above may relate to one embodiment or may relate to several embodiments. The embodiments are not limited to those that solve any or all of the stated problems or those that have any or all of the stated 5] benefits and advantages. It will further be understood that reference to 'an' item refers to one or more of those items.

[00219] The operations of the methods described herein may be carried out in any suitable order, or simultaneously where appropriate. Additionally, individual blocks may be deleted from any of the methods without departing from the scope of the subject matter described herein. Aspects of any of the examples described above may be combined with aspects of any of the other examples described to form further examples without losing the effect sought.

[00220] The term 'comprising' is used herein to mean including the method blocks or elements identified, but that such blocks or elements do not comprise an exclusive list and a method or apparatus may contain additional blocks or elements.

[00221] The term substep' is used herein to refer to a proper substep such that a substep of a step does not comprise all the elements of the step (i.e. at least one of the elements of the step is missing from the substep).

[00222] It will be understood that the above description is given by way of example only and that various modifications may be made by those skilled in the art. The above specification, examples and data provide a complete description of the structure and use of exemplary embodiments. Although various embodiments have been described above with a certain degree of particularity, or with reference to one or more individual embodiments, those skilled in the art could make numerous alterations to the disclosed embodiments without departing from the scope of this specification.

Claims

CLAIMS1. A method for controlling a chemical plant, the method comprising, at a first time interval: simulating the chemical plant using a digital twin, the digital twin comprising: a plurality of connected blocks, each block representing part of a process carried out by the automated chemical plant; for each block, mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block; for each block, one or more second equations describing the part of the process in the block and formed using one or more of first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships, and wherein the digital twin is configured such that, for each block, outputs of the equations are constrained to agree with the mass and energy conservation equations; and triggering control of the chemical plant using results of the simulation; and wherein the method is repeated for subsequent time intervals 2. The method of claim wherein at least one of the blocks has a representation of a thermal boundary to the ambient environment.3. The method of any preceding claim wherein the phenomenological relationships are relationships that take into account one or more of: physical attributes of the chemical plant, chemical attributes of processes in the chemical plant.4. The method of any preceding claim wherein the second equations are formed from first principle relationships and applied by lumping equation parameters together as tunable constants.5. The method of any preceding claim wherein the empirical relationships mathematically relate inputs and outputs of a rate process and are developed and adjusted periodically using error regression comparison of predictions from the digital twin and corresponding measurements from the chemical plant.6. The method of any preceding claim wherein the machine learning derived relationships are adjusted periodically by retraining a machine learning system using measurements from the chemical plant labeled with ground truth values of lumped parameters of the machine learning derived relationships.7. The method of any preceding claim wherein the plurality of blocks comprise at least one block that comprises a representation of one or more of the following processes: a stoichiometric chemical reaction, a dissolution from solid to liquid phase, a precipitation from liquid to solid phase, an absorption from gas to liquid phase, a vaporization from liquid to gas phase, separation of material by phase of matter separation by centripetal forces, separation by filtration.8. The method of any preceding claim wherein the process carried out by the chemical plant is any of production of sodium bicarbonate, production of ammonia, production of metformin, desulfurization of flue gas, carbon dioxide removal from a gas using an amine solvent, the production of hydrogen using solar power.9. The method of any preceding claim comprising using an output of the simulation to change at least one controllable input to the chemical plant in order to compensate for uncontrolled variations in at least one input to the chemical plant 10. The method of any preceding claim wherein, in response to the simulation predicting that a quality of output of the chemical plant can be maintained, changing the at least one controllable input; and in response to the simulation predicting that a quality of output of the chemical plant cannot be maintained, transitioning the chemical plant to a stable state; and in response to the simulation predicting continued deterioration of the uncontrolled variations, moving the chemical plant to a stable state.11. The method of any preceding claim wherein the digital twin comprises, for at least one block, a particle population conservation equation expressed along a particle size coordinate, the particle population conservation equation allowing calculation of changes in particle size distribution as state of the block changes over time and wherein outputs of the second equations are constrained to agree with the particle population conservation equation.12. The method of claim 1 wherein at least one of the blocks of the digital twin represents process in a reaction vessel in the chemical plant, and wherein the mass and energy conservation equations for the at least one block are: a gas subblock mass balance conservation equation which is expressed in words as the rate of change of mass in a gas subblock of the block over time is equal to the sum over species i involved in the process, of the rate of inflow of mass of species (i) minus the sum over species i of the rate of outflow of mass of species (i) plus the sum of the rate of change of mass of species i due to subblock physical and chemical reactions; a bubble phase mass balance conservation equation expressed in the same way as for the gas subblock mass balance conservation equation; a particle population balance conservation equation which is expressed in words as a particle population density distribution partial derivative with respect to time in the block or subblock plus a partial population density derivative of a product of the particle population density distribution and a growth rate, with respect to a particle size coordinate equals a partial derivative of the particle population density distribution with time due to inflow of solid particles minus a partial derivative of the particle population density distribution with time due to outflow of solid particles plus a contribution to the particle population density distribution due to formation rate of new particles via nucleation plus an aggregation kernel that takes allows for the particle population density distribution changes such that a net mass change rate due to agglomeration remains zero over time plus a breakage kernel that takes care of partial particle population density distribution changes such that a net mass change rate due to breakage remains zero over time; a column block energy balance conservation equation which is expressed in words as a product of the sum of the product of each species i mass and specific heat capacity at the block or subblock temperature and the complete derivative of the block or subblock temperature with respect to time minus the change of heat exchange of the block or subblock with its environment equals the total enthalpy change due to the sum of the block or subblock inflow streams minus the total enthalpy change due to the sum of the block or subblock outflow streams plus the sum of the product of each species i mass change and their specific enthalpy at the block or subblock temperature, and a slurry phase mass balance conservation equation expressed in the same way as for the gas subblock mass balance conservation equation.13. The method of any preceding claim wherein the digital twin comprises equations describing heat transfer through walls of pipes and vessels in the chemical plant, and wherein the digital twin is configured to predict the temperature of a primary reaction zone in the chemical plant from weather conditions and measurements obtained from a physical temperature probe in the chemical plant and outside the primary reaction zone.14. The method of claim 13 comprising using the predicted temperature of the primary reaction zone to control an external heat sink on an input stream to the chemical plant and/or to control a thermal jacket on a vessel in the chemical plant.15. The method of any preceding claim comprising detecting downtime of equipment in the chemical plant, selecting a response algorithm from a store of response algorithms according to characteristics of the detected downtime of equipment, and applying the selected response algorithm to control the chemical plant in combination with predictions from the simulation 16. The method of any preceding claim comprising receiving a target for production of an output by the chemical plant and, based on a difference between a predicted output computed by the simulation with the target, controlling the chemical plant so as to reduce the difference 17. The method of claim 16 wherein the target is any one or more of: a quantity of the output, a purity grade of the output, a grade of particle size distribution of the output.18. The method of any preceding claim comprising dynamically comparing a rate measured in the chemical plant with a simulated rate from the simulator, and in response to the comparison meeting criteria, detecting an error in a flow meter in the chemical plant.19. The method of any preceding claim comprising comparing a measured temperature of a process stream in the chemical plant after a heat exchanger to a temperature simulated by the digital twin; and detecting scaling on surfaces of the heat exchanger based on the comparison The method of claim 6 comprising monitoring change in error regression correction factors in order to detect equipment failing faster than other equipment in the chemical plant.21 The method of any preceding claim wherein the digital twin comprises three connected blocks being a first block representing a sodium carbonate and sodium bicarbonate reaction chamber process, a second block representing a slurry recycle, heat sink and sodium bicarbonate production offtake process, and a third block representing a gas recycle and makeup CO2 gas process.22. The method of any preceding claim 21 wherein at least one of the blocks is a chemical reaction block comprising two or more substeps that enable discrete speciation of chemical components, each substep represented using an equation, whereby the equations are solved simultaneously in a given time interval using the digital twin.23. A digital twin for use in controlling a chemical plant, the digital twin comprising a plurality of connected blocks, each block representing part of a process carried out by the chemical plant; for each block, mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block, for each block, one or more equations describing the part of the process in the block and formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships; and wherein the digital twin is configured such that, for each block, outputs of the equations are constrained to agree with the mass and energy conservation equations; and wherein the digital twin has an output for sending simulation results to a SCADA for controlling the chemical plant using results of the simulation 24. A method for designing a chemical plant, the method comprising, at a first time interval: simulating the chemical plant using a digital twin, the digital twin comprising: a plurality of connected blocks, each block representing part of a process carried out by the automated chemical plant, for each block, mass and energy conservation equations that mathematically describe how matter and energy are conserved in the block; for each block, one or more second equations describing the part of the process in the block and formed using one or more of: first principle relationships, phenomenological relationships, empirical relationships, machine learning derived relationships; and wherein the digital twin is configured such that, for each block, outputs of the equations are constrained to agree with the mass and energy conservation equations; and and wherein the method is repeated for subsequent time intervals; and using the digital twin to predict performance of any one or more of: reactor vessels of different dimensions, gas sparging apparatus, insulation, physical characteristics of the chemical plant.25. The method of claim 24 comprising using the digital twin to predict performance of the chemical plant for one or more of: a range of equipment performance parameters of equipment in the chemical plant, a range of weather conditions.