At first, a model was defined that describes the German passenger car population for the next ten years. Therefore, 2019 is considered as a baseline with real data. From 2020–2030, we take recourse to prognosis data, as shown in Section 3.3.

To obtain a sufficient dataset for the 2019 baseline, the vehicle data from the 2017 "Mobility in Germany" (Mobilität in Deutschland, MID) study by the German Federal Ministry of Transport and Digital Infrastructure [17] were categorized into segments and drive train technologies according to Section 3.1. Given the negligible change of 0.2 years in the average vehicle age between 2017 and 2019 [18], we assumed that the overall age distribution of the German vehicle stock did not change significantly. Therefore, a linear transformation of the 2017 data was made to obtain the age distribution of the segments in 2019. Figure 2 presents the resulting age distribution of the German passenger car fleet.

The same data are also available for the vehicle drive train technology and are linked with the year of manufacture and the segment category.

Combined with the number of the total vehicle stock [19], the baseline for 2019 was built.

**Figure 2.** Age distribution of the German vehicle fleet in 2019 (vehicles in scope: 82% of overall vehicle stock).

#### *3.3. German Car Population Development until 2030*

In order to predict the car population until 2030, estimations for new registrations and deregistrations must be made.

#### 3.3.1. New Passenger Car Registrations

To determine new passenger car registrations, the total number and their distribution in the categories defined in Section 3.1 are needed. The annual number of new registrations and their classification regarding the drive train technology were obtained from the "proKlima" scenario of Agora 2018 [20] (Figure 3).

Additionally, an estimation concerning the segments was created based on the trend of the new registration numbers from previous years. For the large van, mini van, upper class, and upper middle class segments, we assumed that the share of new registrations remains constant. In addition, we assumed that the new registration share of segments not considered in this analysis will stay constant, which results in a non-varying share of the overall scope of 82%. A noticeable growth of the sport utility vehicle (SUV) segment has been observed over the last years [21]. The vehicle segment development is visualized in Figure 4.

This trend is expected to continue and is anticipated to reach an annual new registration share of 40% in 2030. Since this value is based on an assumption, a sensitivity analysis is performed in Section 4.4.2.

The growth in the SUV segment was subtracted proportionally from the other segments depending on their share in the reference year 2019. The shares of the large van and mini class segments were not modified, as customers are unlikely to change to an SUV due to the completely different usage profiles.

**Figure 3.** Linear interpolation of the "proKlima" drive drain technology scenario.

**Figure 4.** Segment scenario.

#### 3.3.2. Vehicle Disposal

In this analysis, we assumed that every vehicle that is deregistered from the German transport system, as defined in Section 3.1, will immediately be disposed of, meaning scrapped or recycled. As usage profiles of exported vehicles are unknown, the vehicles are disposed of simultaneously with their deregistration in this model.

For the last four years, a consistent age-dependent deregistration trend can be seen in Figure 5. Therefore, the future deregistration rates are assumed to be identical to those in 2019. The peak after three years can be explained by the end of many leasing contracts, which can contribute to exportation of vehicles. Additionally, the first legally required technical inspection ("Hauptuntersuchung", HU) of the vehicles and the associated phasing out of early defective vehicles take place in this period. This might add to the observed peak. The peaks that occur every two years, especially for older vehicles, are also due to the interval of the mandatory technical inspections in Germany. As technical defects are detected and repairs may become necessary, it is often cheaper to scrap or export the vehicles than to continue their operation.

As explained in Section 3.4, vehicle production is modeled separately from deregistrations, and only for new vehicles. The negligible number of new registrations (negative deregistration) for very old vehicles (probably as vintage vehicles) is ignored.

**Figure 5.** Age-dependent deregistration data for 2016–2019 (based on [22]).

#### *3.4. Modeling Car Population Development*

The development of the car population for different scenarios is calculated using an iterative deregistration and production process. In the first step, the vehicle population (starting with the 2019 population from [22]) for each year of manufacture is reduced according to the deregistration rate for the respective age. Subsequently, the age groups are shifted back one year, and the vehicle production of the considered year is set to maintain a constant total vehicle population. The vehicle deregistration rate is scaled from the real data shown in Figure 5 in order to match the prognosis of newly registered cars between 2020 and 2030 and of the "proKlima" [20] scenario.

In the next step, the modified scenarios (shown in Figure 6) are derived from the baseline by adapting the deregistrations only. First, the deregistration rate is decreased to reflect the reduced registrations due to owners' financial uncertainty during the COVID-19 pandemic. The total vehicle registrations in 2020 are decreased by 23% compared to 2019 [1]. To simulate a subsidy in 2021, the deregistration of vehicles older than 10 years is increased to match the total number of subsidized vehicles. For the baseline scenario, this number is set to the predicted new registration decrease due to the COVID-19 pandemic (0.72 million). In Section 4.4, we discuss the impact of the number of subsidized cars; therefore, additional analyses of twice and four times the baseline amount (1.43 and 2.87 million) were made.

This approach considers the reduced vehicle production in subsequent years after the subsidy is introduced (pull-out effect) and the increased total vehicle production due to early deregisterations compared to a scenario without a subsidy.

**Figure 6.** Deregistration rate distributions of different scenarios.

Figure 7 shows the annual and the accumulated vehicle registrations for different scenarios. The economic uncertainty due to the COVID-19 pandemic reduces the total number of vehicle registrations by forcing owners to keep their vehicles longer. On the contrary, the subsidy leads to earlier replacements and, therefore, increases the new registrations. As already explained for Figure 5, there is a first peak of vehicle decommissioning after three years.

This can also be seen in the new registrations, as the additional vehicles produced due to the subsidy also lead to an increased number of vehicles being taken out of service three years after the subsidy. In turn, this leads to an increase in new purchases that can be seen four years after the subsidy. This time delay is due to the fact that we model the transition from one year to the next, and thus, the ages of the vehicles in the previous year are taken into account for the scrapping and new purchases in the year under consideration.

**Figure 7.** Annual and accumulated vehicle registrations.

#### *3.5. Life-Cycle Emissions*

To perform a cradle-to-grave life-cycle analysis, the production, use, and EoL phases are considered. The focus is primarily on the calculation of greenhouse gas emissions (CO2 equivalent emissions, 20-year global warming potential (GWP 20)). In this analysis, the direct CO2 emissions of the combustion engines in the use phase are set identically to their CO2 equivalent emissions.

The amounts of production and EoL emissions are assumed to be constant for the upcoming years. They are only dependent on the segment and power train type of the vehicle and do not change with the production and disposal date.

In order to analyze the life-cycle emissions of the drive train technologies per segment, specific input data (weight, direct CO2 emissions, fuel, and electric energy consumption) are needed. Complete data for gasoline and diesel vehicles are available. For PHEVs and BEVs, the required input is derived from defined reference vehicles. Currently, in some segments, there are no vehicles available. In these cases, input data are interpolated from adjacent segments.

#### 3.5.1. Production

All production analyses were performed with the Ecoinvent 3.5 database [23]. The cutoff allocation was used and all location settings were set to global.

#### Petrol and Diesel

To perform a life-cycle analysis for gasoline and diesel vehicles, the curb weight was used as the input value (Table 2).

**Table 2.** Weight distribution of vehicles with combustion engines [24].

It was assumed that there is no curb weight difference between gasoline and diesel vehicles. The emissions were calculated with the processes *passenger car production, gasoline* and *passenger car production, diesel* in Ecoinvent. To consider production only, the Ecoinvent data (flow) *manual dismantling of a used passenger car with internal combustion engine* were excluded.

#### Battery Electric Vehicles (BEVs)

In order to perform the production analysis for the BEVs, the curb weight is divided into the battery weight and the remaining vehicle weight. The battery weight is calculated based on the battery capacity. Therefore, a constant energy density for all BEVs is defined.

A representative energy density at a battery packaging level is defined by averaging the values of the VW E-golf (113 Wh/kg) [25] and the 2017 Tesla Model S P100D (160 Wh/kg) [26]. This results in an energy density of 136 Wh/kg, which is used to calculate the battery weight of all BEVs based on the battery capacity. For each segment, a reference BEV is defined based on the number of newly registered vehicles in 2019 [27] (Table 3). In the SUV segment, the Hyundai Kona Electro and the Audi E-Tron both had nearly the same registration numbers in 2019. Since the share of compact SUVs is expected to continue growing [27], the Hyundai Kona Electro is chosen as the reference model.

**Table 3.** Weight distribution of battery electric vehicles (BEVs).

The emissions emitted from the BEV battery pack production were determined with the process *battery production, Li-ion, rechargeable, prismatic*. This process also considers the transportation of the single cells from Peking to Amsterdam by ship and a 1000 km transportation route within Europe by truck.

For the remaining electric vehicles, the process *passenger car production, electric, without battery* [34] was used with the remaining curb weight from Table 3.

#### Plug-in Hybrid Electric Vehicles (PHEVs)

In this analysis, the PHEV is a combination of a conventional and an electric vehicle. Therefore, the PHEV consists of a battery, electric power train components, and the remaining vehicle, which contains all components of the combustion power train. This analysis is based on manufacturers' data from reference vehicles (Table 4).

**Table 4.** Manufacturers' data for plug-in hybrid electric vehicles (PHEVs).

Similarly to the BEV batteries, the PHEV battery weight is calculated based on the battery capacity with a constant energy density for all types of plug-in vehicles. We consider all batteries at the packaging level. Since the PHEV batteries are significantly smaller than the BEV ones, an assumption of the same energy density is not suitable.

Therefore, the mean value of the energy densities of the Kia Optima (75 Wh/kg) [37] and the Mercedes B 250 e (104 Wh/kg) [42] of 89.5 Wh/kg are assumed for all PHEVs.

In addition, the electric power train components, such as an electric motor, a charger, and cables, are considered. Due to lacking data on the actual weight of these components, the weight is scaled with the maximum power of the electric power train. In this analysis, the baseline is made according to the default value in Ecoinvent of a 100 kW power train with a weight of 70 kg. The resulting parameters (Table 5) are used in Ecoinvent for the production analysis.

**Table 5.** Calculated weight distribution of PHEVs.

For the battery pack production of the PHEVs, the same method as that used for the BEVs was applied. The process *electric motor production, vehicle* was used to determine the emitted emissions to produce the components of the electric power train.

The remaining part of the vehicle is assumed to be a gasoline engine car. Therefore, the process *passenger car production, petrol* was applied to analyze the production.

#### 3.5.2. Use Phase

The emissions emitted during the use phase are composed of the emissions based on fuel and electric energy consumption, the well-to-tank emissions, and the emissions due to maintenance of the vehicle and the road infrastructure.

In 2019, the average annual mileage of passenger cars was 13,602 km [43]. In this model, it is assumed that, in the future, every car will have the 2019 average annual mileage.

#### Combustion Engine Vehicles

The use phase emissions of the gasoline- and diesel-powered vehicles are dependent on the vehicle's segment and age. Based on these criteria, the KBA [44] provides a database of the average direct CO2 emissions and fuel consumption, which are displayed in the supplementary materials in Table S1. In this work, they are assumed to be identical to the CO2 equivalent emissions.

#### Battery Electric Vehicles (BEVs)

In this analysis, the BEVs cause indirect CO2 emissions because of their consumed electric energy during the use phase.

For each segment, a reference BEV is defined based on the number of newly registered vehicles in 2019 [27]. If a model in this segment would have a dominant share of the newly registered cars, it is considered as the reference BEV. Otherwise, the two models with the highest number of newly registered vehicles in the segment are selected, and a weighted average of the energy consumption (weight factor) according to their new registration numbers is built (Table 6).

In 2019, there were not any BEVs newly registered in the upper middle class or large van segments. Since the market for electric vehicles will grow [20], we assume that future electric vehicles will also be available in these segments. Based on the curb weight of the combustion engine vehicles, we estimate that upper middle class and large vans have the same energy consumption. The electric energy consumption of upper middle class is assumed as an average of those of the middle class and upper class.

#### Plug-In Hybrid Electric Vehicles (PHEVs)

To determine comparable CO2 emissions from plug-in hybrids, a combination of the direct exhaust emissions and the energy consumption of the electric drive train is made (Table 7).

The manufacturers' data for the reference vehicles either rely completely on the New European Driving Cycle (NEDC) or were measured with the Worldwide Harmonized Light Vehicles Test Procedure (WLTP) and were transferred back to the NEDC. The NEDC considers half of the driving cycle using the combustion engine and the other half using the electric power train. In these data, the electric energy consumption is considered emission-free. In order to correct this effect, we add the appropriate indirect emissions for electric power consumption.

**Table 7.** CO2 emissions and energy consumption of PHEVs.

#### Indirect Emissions

In this model, the consumed electric energy causes indirect emissions based on the specific carbon dioxide emissions of the German electricity mix. In 2020, the generation of one kWh electric energy is predicted to cause 432 g CO2 eq. emissions [48].

Pehnt et al. [48] calculated 318 g CO2 eq. emissions for the electricity mix in 2030; the value here decreases linearly for the years 2020 to 2030. The well-to-tank factors are used to determine the upstream emissions of consuming fuel. Here, Schallaboeck et al. [49] stated that 685 g or 408 g CO2 equivalent emissions produce one liter of gasoline or diesel, respectively. The data for the fuel consumption are displayed in the supplementary materials in Table S2. All PHEVs are assumed to have a gasoline engine.

#### Maintenance

During the use phase, the replacement of spare vehicle parts and the maintenance of the road network are considered. To calculate the vehicle maintenance, the flow *passenger car maintenance* was used in Ecoinvent and was adapted according to the vehicle's weight and power train type [50]. We followed the same approach as in the Agora (2019b) [51] "basic scenario", assumed for each vehicle a lifespan of 150,000 km, and did not consider any battery exchange for BEVs and PHEVs.

In order to determine the maintenance of the road network, the flow *road maintenance* is used. This introduces a weight-dependent contribution during the use phase.

#### 3.5.3. End-of-Life

In this model, it is assumed that every deregistered vehicle is eventually disposed of. The equivalent CO2 emissions due to this process were determined in Ecoinvent using the flows *manual dismantling of used passenger car with internal combustion engine*, *treatment of used glider, passenger car, shredding*, *treatment of used internal combustion engine, shredding*, *treatment of used powertrain for electric passenger car, dismantling*, and *market for used Li-ion battery* according to [34].

#### *3.6. Subsidy Concepts*

Three subsidy concepts are introduced in this chapter. They reflect aspects of the current discussions about the possible subsidy scope. Höpfner et al. [3] showed that 84% of the subsidy from the "environmental bonus" in 2009 was used to purchase small cars (mini class, small class, and compact class). Therefore, 84% of all new registrations due to the subsidy will happen in the mini class, small class, and compact class categories. For this distribution, the predicted new car registration shares for 2021 are used. All segment shares are evenly scaled in order to reach an 84% share in the small cars category. The impact of this effect is discussed in Section 4.4.1.

The distribution shown in Figure 8 is used as a baseline for all subsidy concepts. In Section 4.4.1, the influence of this distribution is discussed.

Similarly to the "environmental bonus" in 2009, a mandatory requirement to earn the subsidy is a minimum vehicle age (in this case, 10 years).

**Figure 8.** Subsidized vehicle distribution of the "broad funding" concept.

#### 3.6.1. "Broad Funding" Concept

In this concept, the subsidy is not restricted to a certain vehicle type. Except for the higher demand for the mini, small, and compact class segments, the same composition of newly registered cars as predicted without the subsidy will occur. The subsidy only increases the total number of newly registered and deregistered cars. For every newly registered car due to the subsidy, a car in the same segment, manufactured in 2010 or earlier, is disposed of. Vehicle purchases in higher segments are not subsidized.

#### 3.6.2. "Innovation" Concept

With this concept, only BEVs are subsidized. For the baseline model, the same distribution of subsidized cars as in Section 3.6 is considered, except only BEVs will be newly registered due to the subsidy.

#### 3.6.3. "Downsizing" Concept

Lastly, a concept is shown that only allows a subsidy if the new car is at least one segment smaller (defined by the segment's CO2 emissions) than the one traded in. The respective segment shifts in the "downsizing" concept are displayed in Table 8.

As per the CO2 emission-based downsizing definition, the SUVs and mini vans are shifted to the middle class segment.

#### **4. Results**

The goal of this work is to analyze the impact of a subsidy on the CO2 eq. emissions of the German passenger car fleet. This is achieved by considering categorized vehicle emissions and fleet composition data.

#### *4.1. Life-Cycle Analysis*

The life cycle of a vehicle consists of production, maintenance, use phase, and EoL processes. For each element, the respective emissions are calculated per segment and per drive train technology.

#### 4.1.1. Production

Using the process described in Section 3.5.1, an analysis of the production emissions for each vehicle is made. The results for each vehicle class and drive train technology are shown in Figure 9.

The production of BEVs and PHEVs results in higher CO2 equivalent emissions than the production of conventional vehicles. This phenomenon emerges in every segment and can be attributed, in addition to the higher vehicle weight of BEVs and PHEVs, to the greater effort of producing batteries and electrical parts. The BEV production effort increases substantially for larger segments as more battery capacity is implemented. This effect is less pronounced for PHEVs.

#### 4.1.2. Maintenance

The emissions from the vehicle and road maintenance per driven kilometer (Figure 10) scale mainly with the vehicle weight.

**Figure 9.** Production emissions per vehicle, segmented by class and drive train technology.

**Figure 10.** Maintenance emissions.

#### 4.1.3. Use Phase

The CO2 equivalent emissions depend on the vehicle specifications. In addition to the direct exhaust emissions, vehicles with combustion engines cause indirect emissions via the procurement of the fuel. Figure 11 shows the combined specific CO2 eq. emissions per kilometer in 2019.

**Figure 11.** Composition of CO2 eq. emissions per driven kilometer.

Vehicles with a gasoline engine have the highest specific emissions across every segment, followed by diesel engine vehicles. The BEVs have the lowest specific emissions of the considered drive train technologies in all segments.

#### 4.1.4. End of Life

Due to the high recycling effort for the batteries, we see increased EoL emissions for the BEVs and PHEVs compared to the conventionally powered vehicles (Figure 12). The higher battery capacity in larger segments leads to higher EoL emissions for the BEVs compared to the PHEVs. Again, the vehicle weight is an important driver of EoL emissions.

**Figure 12.** End-of-life emissions per vehicle, segmented by class and drive train technology.

#### *4.2. Subsidy's Impact on Fleet Emissions for 2019–2030*

Using a combination of the results from Section 4.1 and the German car population scenario from Section 3.3, a prediction of the CO2 equivalent emissions emitted from the German passenger car system until 2030 is made. The emissions decrease in 2020 due to reduced vehicle production (Figure 13). As shown in Figure 7, the vehicle production increases in the subsequent years due to the subsidy. Therefore, higher emissions are expected before the vehicle fleet benefits from the increased efficiency that leads to reduced annual emissions.

The cumulative CO2 eq. emissions of the baseline scenario without a subsidy (Figure 13) are compared to the cumulative emissions of the various subsidy concepts. Each concept is calculated with 0.72 million subsidized vehicles. The developed model enables the analysis of any number of subsidized vehicles.

Figure 14 shows the comparison of the different subsidy concepts relative to the baseline scenario. As a result, the accumulated differences of the CO2 equivalent emissions are shown for the respective scenarios.

The cumulative difference shows significant additional emissions in 2021 for all subsidy concepts. The effect displays the higher number of new vehicles and the related production emissions for this year. The emission savings are dependent on the subsidy concept.

The "innovation" subsidy clearly lowers the CO2 equivalent emissions of the German passenger car transport system until 2030 the most out of all considered concepts. Even with the highest CO2 equivalent emissions in 2021, the "innovation" concept reaches the break-even earliest in the year 2025. The additional CO2 equivalent emissions of the "innovation" concept in 2021 are 8.6% higher than the "broad funding" and 20.8% higher than the "downsizing" concepts. This effect can be attributed to the higher production effort for electric vehicles, as shown in Figure 9.

**Figure 13.** Life-cycle CO2 eq. emissions of the baseline scenario without a subsidy.

**Figure 14.** Cumulative difference of CO2 eq. emissions compared to the scenario without a subsidy.

Despite having the highest CO2 equivalent emissions in 2021, the accumulated CO2 equivalent emission savings until 2030 of the "innovation" concept are 285.7% higher than the "downsizing" and 576.7% higher than the "broad funding" concepts. Compared with the scenario without a subsidy, the "innovation" concept reaches a cumulative difference of −7.56 × 106 t CO2 equivalent emissions until 2030.

In the "downsizing" concept, compared to the "broad funding" concept, smaller vehicles are produced. Therefore, in 2021, the CO2 equivalent emissions of the "downsizing" concept are lower than the emissions in the "broad funding" concept. In addition, the "downsizing" concept reaches the break-even point in 2027, one year earlier than the "broad funding" subsidy concept.

#### *4.3. Sensitivities*

Due to the complexity of the model, dependencies on key parameters are reviewed and their respective impacts on the subsidy concepts are determined.

#### *4.4. Number of Subsidized Cars*

As shown in Section 3.4, we assume that the number of subsidized cars matches the predicted registration decrease in 2020 due to the COVID-19 pandemic. To determine the impact of the number of subsidized cars (∼0.7 million), this value is doubled (∼1.4 million) and quadrupled (∼2.8 million).

As expected, the cumulative difference of CO2 eq. emissions in 2030 is proportional to the number of subsidized cars, which can be seen in Figure 15.

**Figure 15.** Number of subsidized cars.

#### 4.4.1. Vehicle Distribution

As discussed in Section 3.6, the "environmental bonus" in 2009 was primarily used to purchase smaller vehicles. This can be attributed to the fairly high share of the subsidy compared to the purchase price.

It would be conceivable to design the subsidy in a way that the amount is determined by a fixed percentage of the purchase price. It is assumed that this leads to a vehicle distribution of the subsidized cars proportional to the predicted new passenger car registrations in 2021.

Figure 16 shows the difference between a scenario proportional to 2021 and a scenario shifted to smaller cars. The modified distribution of subsidized vehicles has an impact on the results.

As shown in Table 9, all subsidy concepts are less efficient when using the distribution proportional to 2021. In this context, the "broad funding" concept even increases the CO2 eq. emissions compared to the no-subsidy baseline. As a result, subsidizing smaller cars significantly increases the efficiency of the subsidy.

**Table 9.** Difference of cumulative CO2 eq. emissions compared to the scenario without a subsidy in 2030.

**Figure 16.** Subsidized car distribution for two scenarios under investigation.

#### 4.4.2. New Vehicle Registrations

#### Drive Train Technologies

To predict the German car population until 2030, we rely on Section 3.3 and on the "ProKlima" scenario stated in Agora (2018) [20] to describe the trend of the drive train technologies. To estimate the impact of the subsidy for a baseline scenario with a higher share of newly registered PHEVs and BEVs, an analysis with the "ProKlima Plus" scenario is made.

In Table 10, the values for the drive train technology development of the different scenarios can be seen. To determine the share of newly registered cars for each year, a linear interpolation between the displayed values is made.

**Table 10.** Annual share of newly registered cars.

Figure 17 shows the cumulative CO2 eq. emission difference of the "innovation" concept with the two different assumptions for the annual newly registered cars until 2030. It is shown that the subsidy loses efficiency if the baseline share of PHEVs and BEVs rises. The cumulative CO2 eq. emission difference will be smaller due to the overall more efficient German passenger car fleet with the "proKlima Plus" scenario.

**Figure 17.** Cumulative CO2 eq. emission difference between the "no subsidy" baseline and the "innovation" concept using the "proKlima" and "proKlima Plus" registration predictions.

#### Segments

In Section 3.3, an assumption of the development of the segment composition of newly registered cars is introduced. Here, we assume that the SUV segment will reach a 40% share of the annual newly

registered cars in 2030. Since this parameter is uncertain, analysis with either a 30% or 50% proportion of SUVs of the newly registered cars in 2030 is made. For the "innovation" concept, the 30% SUV share scenario generates 2.40% less CO2 eq. emissions savings compared to the 40% SUV share baseline due to the overall more efficient car fleet. On the contrary, the cumulative CO2 eq. emissions savings of the "innovation" concept in the 50% SUV share scenario are 0.84% higher compared to the baseline.

This shows that the results of the study are robust against this parameter because the trend and the statement of the results are not changed by a variation of the SUV share.

#### **5. Discussion**

The absolute CO2 eq. emissions savings are proportional to the number of subsidized cars. For the results shown, we set the subsidized cars in the model to about 0.72 million cars. This is more conservative than the actual replacements of 1.95 million [3] during the 2009 "environmental bonus".

Our model is able to deal with a variable number of cars, as shown in the sensitivity analysis in Section 4.4. The actual number of subsidized cars would depend on micro- and macroeconomic factors, as well as technology-specific concerns, with different scenarios needing different funding to reach the same total number of replaced vehicles. Additionally, an investigation of the basic conditions, e.g., charging facilities, is needed to match the implementation of the number of subsidized vehicles. However, since we are not conducting an economic analysis, quantifying these factors reliably is out of the scope of this work. Our model can be easily updated to calculate the effects of arbitrary numbers of subsidized vehicles. In case the 2020 program would have the same impact on car sales, the absolute savings grow accordingly.

In the absence of detailed prognosis data for the technological improvements, we assume no efficiency gains regarding the CO2 eq. emissions. The trend of the combustion engines indicates a saturation for the CO2 efficiency [44]. Assuming that the electric drive train technologies will still gain efficiency [52] and more electricity will be generated from renewable sources, the total savings for the "innovation" concept would be higher.

Since only the direct CO2 emissions for the driving during the use phase of combustion engines are available, the actual CO2 eq. emissions are slightly higher. As this effect is not relevant to BEVs, since CO2 eq. emissions are considered for the electricity production, the "innovation" concept would generate marginally higher savings when considering all emissions from the combustion engines.

This paper shows that all considered subsidy concepts eventually have a positive effect on reducing the CO2 equivalent emissions of the German passenger car fleet. This corresponds with both Höpfner et al. [3] and Lenski et al. [13], who stated that the 2009 car subsidy had a positive effect on the environmental performance of the various car fleets.

Contrary to Klößner et al. [8], the long-term effect on the CO2 emissions is positive, as windfall gains (subsidized purchases that would happen without the subsidy) are not included and BEVs and PHEVs are considered. We only expect additional CO2 emissions to occur immediately after the subsidy is introduced due to the increased production, which is overcompensated by savings during the use phase.

To calculate the LCA emissions, we clustered the vehicle fleet into segments to reflect the actual distribution. Agora (2019) [53] used only one reference vehicle for each drive train technology. For comparable segments, the production emissions of the BEVs show similar results. In this work, the gasoline and diesel vehicles result in higher production emissions compared to those in Agora (2019) [53]. This can be attributed to a deviation in parameter settings.

#### **6. Conclusions**

In order to calculate the cumulative CO2 eq. emissions of the German passenger car fleet from 2019–2030, a granular model was developed. Three different subsidy concepts were introduced and the impact on the CO2 eq. emissions was determined.

Subsidizing the German passenger car system due to the COVID-19 pandemic shows long-term CO2 emissions savings for nearly all investigated scenarios. Only in the sensitivity analysis, the "broad funding" concept leads to slightly increased emissions when using a vehicle distribution proportional to that of 2021. Taking the sensitivity analyses into account, the "innovation" concept shows the most significant emissions savings in the German passenger car system.

Considering the time period 2019–2030 and a total number of 0.72 million subsidized vehicles, the "innovation" concept generates about 7.56 million t less CO2 eq. emissions compared to the scenario without a subsidy. This equates to 0.46% of the total CO2 eq. emissions of the addressed segments in this period. The "downsizing" and "broad funding" concepts create savings of 0.16% and 0.08%, respectively.

In 2021, the increased vehicle production leads to higher CO2 emissions for all subsidy scenarios compared to the scenario without a subsidy. The ecological break-even is reached in 2025 for the "innovation" concept, in 2027 for the "downsizing" concept, and in 2028 for the "broad funding" concept.

If, from an economic and political point of view, a subsidy program for passenger cars in Germany is considered to be desirable, we clearly recommend the exclusive funding of BEVs because the "innovation" concept achieves the highest positive climate impact at the earliest time.

#### **7. Outlook**

To further investigate the environmental impact of passenger car subsides in Germany, additional greenhouse gases in the use phase and air pollutants must be considered. It seems reasonable to add greenhouse gases with a high GWP. Due to the current discussion on driving bans in German city centers, the subsidy's impact on nitrogen oxides would enhance the model's result.

It is conceivable to examine further subsidy concepts, such as focusing on PHEVs or particular segments.

In the current model, the technological status of the vehicles is assumed to remain at the level of 2019. The indirect emissions of only BEVs and PHEVs will decrease until 2030 due to the slight reduction of the specific carbon dioxide emissions of the German electricity mix. This means that the vehicle weight and direct exhaust emissions of each newly registered vehicle will remain constant in the current model.

In our model, we define a reference vehicle to describe the BEVs and PHEVs due to a lack of data provided by the KBA. To determine the detailed PHEV and BEV model parameters more precisely, the data of all vehicles available need to be consolidated.

As Höpfner et al. [3] notes, the vehicle subsidy in 2009 led to reduced used car sales. This was not implemented in the current model and might affect the disposal rate and the number of newly registered cars.

This paper did not investigate the economic impacts of subsidy scenarios. A detailed market analysis is needed to estimate consumers' buying behavior. To improve the model, it is necessary to further examine the subsidy's impact on the newly registered cars and to address economic effects, like on-top sales and windfall gains.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2071-1050/12/23/ 10037/s1: Table S1: Direct CO2 emissions per kilometer, Table S2: Fuel consumption per 100 km.

**Author Contributions:** Conceptualization, L.H. and A.G.; methodology, M.S. and L.H.; software, A.M.S., L.H., and M.S.; validation, A.M.S., A.G., M.S., L.H., and D.G.; investigation, M.S., A.G., A.M.S., and L.H.; resources, A.M.S. and D.G.; data curation, M.S.; writing—original draft preparation, M.S.; writing—review and editing, A.M.S., A.G., L.H., and M.S.; visualization, L.H. and M.S.; supervision, D.G.; funding acquisition, D.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), grant number 398051144, project title: "Analysis of strategies to fully de-carbonize urban transport".

**Acknowledgments:** We acknowledge the support from the German Research Foundation and the Open Access Publication Fund of TU Berlin.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Removal of Basic Brown 16 from Aqueous Solution Using Durian Shell Adsorbent, Optimisation and Techno-Economic Analysis**

**Yashni Gopalakrishnan 1, Adel Al-Gheethi 1,\*, Marlinda Abdul Malek 2,\*, Mawar Marisa Azlan 1, Mohammed Al-Sahari 1, Radin Maya Saphira Radin Mohamed 1,\*, Sadeq Alkhadher <sup>1</sup> and Efaq Noman 3,4**

Received: 2 September 2020; Accepted: 27 September 2020; Published: 27 October 2020

**Abstract:** Azo dyes including C. I. Basic Brown 16 (BB16) are one of the coloured organic compounds that have adverse effects on human health and the environment. The current work aims to optimise the adsorption of C.I BB16 in aqueous solution using durian (Durio zibethinus murray) shell as a low-cost green adsorbent. Durian shell was characterised by Fourier transform infrared spectroscopy (FTIR) and scanning electron microscopy (SEM). The adsorption process was optimised with response surface methodology (RSM) based on pH (4–8), time (30–240 min), durian shell dosage (0.1–1.0 g/L) and initial concentration of C.I BB16 (10–20 ppm). The removal efficiency was determined based on the reduction of chemical oxygen demand (COD) and the decolourisation of C.I BB16. The techno-economic analysis was described in the current work to know the economic feasibility of durian shells as an adsorbent. The SEM images showed that durian shell adsorbent has a smooth surface with no pores. FTIR spectra confirmed the presence of -C-O, =C–H, C=C, -C-O-C and O-H bonds in durian shell. Maximum decolourisation (77.6%) and COD removal (80.6%) for C.I BB16 was achieved with the interaction between pH, time and adsorbent dose and initial concentration of C.I BB16. The optimal operating factors for adsorption of C.I BB16 recorded at pH 8, time (30 min), durian shell dosage (1 g/L) and 15 mg /L of C.I BB16 concentrations were 77.61 vs. 74.26 (%) of C.I BB16 removal and 80.60 vs. 78.72 (%) of COD removal with an R<sup>2</sup> coefficient of 0.94 at *p* < 0.05. The specific cost of durian shell coagulant production is USD 172.71 per ton which is lower than the market price of honeydew peels-activated carbon (HDP-AC) (USD 261.81) and the commercial market price of activated carbon which is USD 1000.00/tons. These findings indicated that the durian adsorbent provides alternative methods for treating hair dye wastewater. These findings indicated that durian shells have a high potential for the adsorption of C.I BB16 in aqueous solution.

**Keywords:** Basic Brown 16; adsorption; durian shell; mechanism; techno-economic analysis

#### **1. Introduction**

Azo dyes are coloured organic compounds bearing the functional group R−N=N−R , in which R and R are aryl or heteroaryl. More than 0.7 million tons of these compounds are produced every year from which 10–15% enters into the environment through industrial wastewater, especially textile industries [1]. The improper disposal of these wastewaters is the main problem of environmental pollution because of its high toxicity and stability to withstand degradation even under extreme circumstances [2]. The main concerns related to the presence of azo dyes in water ecosystems are its effect on water transparency, photosynthetic activity reduction and the harmful consequences on aquatic flora and fauna. Moreover, most of the azo dyes have been documented to be mutagens or carcinogens to human beings [3]. C. I. BB16 is one of the most applicable dyes in the textile and biological dyeing industries. It contains (-N=N-) as chromophore connected to aromatic systems with lateral groups, including CH3, which can be very harmful and resistant to conventional treatment processes. Although C.I BB16 has caused health hazards because of which it has been prohibited in many areas, it is still widely practised in several countries [4].

Several treatment technologies for removing azo dyes from industrial wastewater have been used such as coagulation-flocculation [5], chemical and electrochemical oxidation [6], ozonation [7], ion exchange [2] and anaerobic-aerobic system [8]. However, most of the azo dyes are still persistent in the treated wastewater because of their high stability. Therefore, the removal of azo dyes required an advanced removal technology such as adsorption which is one of the most favourable processes in comparison to other practices because of its simplicity, practicality, productivity and inexpensiveness [9]. Adsorption is also known as physical attachment or binding of ions and molecules onto the surface of other molecular substances [10]. The exploration for new adsorbents has been concentrated on biomaterials because of their high availability, eco-friendliness, economic feasibility and source i.e., from renewable materials. These resources are locally and easily accessible in considerable quantities in the agroindustry. Thus, they are typically inexpensive and have little or no economic significance. For instance, these resources are wastes with no industrial or agricultural use, with their accumulation being problematic in many cases [11]. The adsorbents comprise proteins, polysaccharides, cellulose and lignin, which contain various functional groups such as phenolic and amino groups that are accountable for pollutant adsorption. Among several plants, pineapple plant stem [12], leaves of Prunus Dulcis [13], tomato waste and apple juice residue [14], grape peel [15], *Thapsia transtagana* stems [16], mushroom [17] have been reported to produce biocompatible adsorbents.

In this study, Durio Zibethinus Murray (durian) shell was utilised for the adsorption of C.I BB16 because of its richness in proteins, polysaccharides, cellulose and lignin that contain various functional groups such as phosphate, carboxyl, sulphate, phenolic and amino groups that are accountable for azo dyes adsorption [18]. Moreover, the durian shell contains high fibre content that contributes to the adsorption process [19]. Durian peels are broadly obtainable as a waste from the durian processing business at a cheaper rate in comparison to other natural supplies. This is because only one-third of durian is digestible (20–30% flesh), whereas the shell and the seed (55–66%) are generally disposed of [20]. The productivity of durian in Malaysia alone was 3800 metric tons in the year 2016 and 22,000 metric tons in the year 2020. With every 10 tons, more than 7 tons of durian shells are produced [21]. Although durian shell as adsorbent was investigated by previous researchers [21,22], no comprehensive study has been conducted on the optimisation of adsorption of C.I BB16 by using a response surface methodology (RSM) and this emphasises the novelty of the current work.

Figure 1 presents the bibliometric analysis of 43 papers published on the durian shell according to the Scopus database using specific keywords such as "durian" AND "adsorption". It was noted that the durian shell was used as an adsorbent for removing methylene blue, basic dyes, brilliant green, basic green four dye, crystal violet, gentian violate, acid green 25 as well as heavy metals such as chromium, lead, zinc from different wastewater sources. The main factors affecting the adsorption process included pH, time, the concentration of pollutants, moisture and particle size of the durian

shell powder. The durian shells have been subjected to acid treatment to prepare charcoal, biochar and carbon active which are the most common adsorbents.

**Figure 1.** Bibliometric analysis for 43 paper published on the durian shell according to the Scopus database using specific keywords included "durian" AND "adsorption".

The factors affecting adsorption need to be optimised to reinforce the process of adsorption. Response surface methodology (RSM) is a tool that uses the mathematical and statistical approach to develop, optimise and enhance products and industrial processes. The main aim of RSM is to lessen the number of experimental runs and attain maximum data output. It allows changing numerous variables at the same time to plan fewer number of experiments needed to optimise the desired process [23]. RSM can be used to study the effects of independent variables, multiple factors and their interaction effects [24]. The effects of two or more variables are determined by factorial design. Meanwhile, central composite design (CCD) with RSM is used for the optimisation process. Besides, the application of CCD prevents deceiving conclusions and grant conclusions that are appropriate for several experimental conditions. Only one study has been optimised for the durian seed prepared as activated carbon for removing C.I methylene blue [25]. However, RSM is not used as a modelling and optimisation tool for the adsorption of C. I BB16. Moreover, techno-economic analysis is used for evaluating the potential of agriculture wastes as an alternative adsorbent for the commercial activated carbon (AC). No study has investigated the final price of the durian shell compared with that of commercial products in the market. The current price of 1 ton AC from coal ranges from USD 1190 to 16343/tons AC [26]. Therefore, the current study aims to optimise the factors affecting the adsorption of C.I BB16 by durian shell as an adsorbent by RSM. Moreover, techno-economic analysis for the preparation and application of durian shell adsorbent was conducted and compared with that of commercial AC products.

#### **2. Materials and Methods**

#### *2.1. Preparation of Durian Shell as an Adsorbent*

The durian shells were gathered locally from the neighbourhood shops. The durian shell was rinsed two to three times by water for the removal of dirt and impurities followed by open-air drying for 48 h. Thereafter, it was oven-dried at 150 ◦C for 24 h. The dried durian shell was then ground and sieved to the particle size <150 μm. The ground durian shell was kept in the airtight plastic container for the adsorption studies [22]. It was noted that 1 kg of durian shell generates 500 g of the powder adsorbent with particle size <150 μm.

#### *2.2. Characterisation of Adsorbent*

#### 2.2.1. Field Emission Scanning Electron Microscope (FESEM)

The surface texture and morphology of the durian shell was analysed by field emission scanning electron microscopy (FESEM; Ion Company, New York, NY, USA). The sample was prepared by mixing two drops of 100% ethanol with 0.1 g durian shell on clean silicon wafers and then dried at 50 ◦C for 2 min to eliminate the solvent, coated with gold using the sputtering technique, fractured and observed at a beam energy of 20 kV.

#### 2.2.2. Fourier Transform Infrared (FTIR)

The functional groups in the durian shell were identified by Fourier transform infrared (FTIR) spectroscopy using Nicolet iS50 (Thermo Scientific, Waltham, MA, USA) with a scanning range of 400–4000 cm−<sup>1</sup> using attenuated total reflection (ATR) technique. First, the crystal plate (sample holder) was cleaned with acetone. The durian shell sample (5 mg) was placed directly on the crystal plate and then pressed using the metallic pressing device and scanned to get the IR spectrum [27].

#### *2.3. Preparation of C.I BB16 in Aqueous Solution*

C.I BB16 was supplied by Sigma-Aldrich (M) Sdn Bhd, Malaysia, and used as received without further purification. A stock solution of 1000 mg/L was prepared by dissolving the appropriate amount (1000 mg) of C.I BB16 in a litre of deionised water. The working solution was prepared by diluting the stock solution with deionised water to give the appropriate concentration of the working solutions [28].

#### *2.4. Experimental Set-up*

The experimental set-up in this study includes the preparation of durian shell and C.I BB16 solution, independent factors and optimisation process. Factorial complete randomised design (CRD) (3 × 3 × 3) in triplicate was executed to investigate the independent factors that affect the adsorption of C.I BB16. Four independent factors that are pH (4–8), time (30–240 min), durian shell dosage (0.1–1.0 g) and initial C.I BB16 concentration (10–20 ppm) were investigated. The range for each independent factor was chosen based on previous investigations [29,30].

#### 2.4.1. Adsorption of C.I BB16

The adsorption of C.I BB16 experiments were performed using a 250 mL Erlenmeyer conical flask that holds 100 mL of an aqueous solution of C.I BB16. To ascertain the pH condition, the pH was modified deliberately for the solution, ranging between 4 and 8. These were conducted by using 0.1 M sodium hydroxide (NaOH) and 0.1 hydrochloric acid (HCI) to acquire the desired pH value [31]. The mixture was stirred using the orbital shaker at 150 rpm and the contact time was specified between 30 min and 4 h. Treated C.I BB16 samples were filtered using the membrane filter (0.2 μm) [16]. The efficiency of the durian shell in the adsorption of C.I BB16 was investigated by determining the absorbance of C.I BB16 sample before and after each experimental run using UV-Vis spectrophotometer

(DR6000) at 661 and 620 nm. The adsorption of C.I BB16 efficiency was determined by calculating the decolourisation according to Equation (1) [32]

$$\begin{array}{rcl}\text{Decolorisation (\%)} &=& \frac{A\_0 - A}{A\_0} \end{array} \tag{1}$$

where *A*<sup>0</sup> is the initial absorbance of C.I BB16 in samples before the adsorption and *A* is the absorbance of C.I BB16 after the adsorption—moreover, the mineralisation of C.I BB16 was confirmed by COD analysis, according to Gonçalves et al., [11]. *COD* reduction was calculated using Equation (2)

$$\text{COD (\%)} = \frac{\mathcal{C}\_0 - \mathcal{C}\_f}{\mathcal{C}\_0} \tag{2}$$

where *C*<sup>0</sup> is the initial *COD* concentration (mg/L) and *Cf* is the final *COD* concentration (mg/L) after the adsorption.

#### 2.4.2. Optimisation of Adsorption of C.I BB16

The experimental set-up of the current work consists of C.I BB16 solutions and durian shell adsorbent as well as pseudo-first and second-order kinetic model for adsorption process using response surface methodology (RSM). The factorial complete randomised design (CRD) (4 × 2) in triplicates was used to study the optimal factors affecting the adsorption process with four (4) independent factors ((pH 4–8), *x*<sup>2</sup> (time, 30–240 min), *x*<sup>3</sup> (adsorbent dosage (0.1–1 g/L), *x*<sup>4</sup> (initial concentration of C.I BB16) (10–20 mg/L) and two dependent variables *y*<sup>1</sup> (decolourisation of C.I BB16) (%), *y*<sup>2</sup> (COD removal) (%). The independent factors were selected and used for designing the experimental runs using Design Expert 6.0.10, central composite design (CCD) (Stat-Ease Inc., Minneapolis, MN, USA) software that suggested 25 experimental runs. The experimental runs (one run recommended) were carried out to authenticate the efficiency of the results of the evaluated elements. Four independent variables including pH, time, durian shell dosage and initial concentration of C.I BB16 were optimised for adsorption of C.I BB16. The optimisation experiments were also operated in three replicates to enhance the estimated data accuracy by inferring that no bias and systematic error exist. Moreover, the variables' ranges (high and low) in this study are depicted as (low, medium and high, being coded as −1, 0 and +1). The efficiency of adsorption of C.I BB16 was verified by the percentage of the removal (*y*) which was denoted as the responses of the parameters of the experimental design. The parameters examined for adsorption of C.I BB16 were decolourisation (%) (*y*1) and COD removal (%) (*y*2) and the working variables were *x*1, *x*2, *x*<sup>3</sup> and *x*4. The quadratic model used to get the optimisation of the responses as a function of the independent process variables is as follows:

$$\|y\| = \beta 0 + \sum\_{i=1}^{k} \beta\_i \mathbf{x}\_i + \sum\_{i=1}^{k} \beta\_{ii} \mathbf{x}\_i^2 + \sum\_{i=1}^{n} \sum\_{i$$

where *y* is the expected response for the adsorption, while β0, β*i*, β*ii* and β*ij* are the linear model coefficients, *xi* is the coded independent factors, and *xi* is the independent factors in a coded form, while *k* represents independent variable numbers.

The relationship between the experimental data and the coded variable is represented as:

$$\varepsilon\_{i} = \varepsilon\_{i} \frac{[\text{HL} + \text{LL}]/2}{[\text{HL} - \text{LL}]/2} \tag{4}$$

where *xi* represents the coded variable, ε*<sup>i</sup>* is the experimental laboratory data, HL and LL represent the maximum and minimum values of the independent factors respectively.

#### *2.5. Techno-Economic Analysis for Hair Dye Greywater Treatment Using Durian Shell Adsorbent*

The annual time for the fill operation including preparation and application is 5760 h/year (which is equivalent to 240 operating days). The TCI for a proposed preparation and application processes includes WCC and FCE (Equation (5)) [33,34].

$$TCI = \text{WCC} + \text{FCE} \tag{5}$$

where *TCI* is the total capital investment, *WCC* is the working capital cost and *FCE* is the fixed capital estimation.

The annual operation cost (*AOC*) for the preparation of durian shell as adsorbent and then the application of the adsorbent in the hair dye greywater and the costs of raw materials (*CRM*), waste generated from the preparation process (*CWG*), utilities (*CU*) and extra cost (*CE*) were calculated on an annual basis according to Equation (6).

$$AOC = \mathcal{C}\_{RM} + \mathcal{C}\_{WG} + \mathcal{C}\_{U} + \mathcal{C}\_{E} \,\mathrm{(}\,\mathrm{s}\tag{6}$$

According to the mathematical model designed by the specific cost of durian shell adsorbent production *CZ* (kg/USD) is calculated based on the annual costs of capital *CC*, utilities cost *CU*, raw material cost *CRM* and extra cost *CE* divided by the annual durian shell adsorbent production (Equation (7)) [35].

$$\mathcal{C}\_{\mathcal{Z}} = \frac{(\mathcal{C}\_{\mathcal{C}} + \mathcal{C}\_{\mathcal{U}} + \mathcal{C}\_{\mathcal{RM}} + \mathcal{C}\_{\mathcal{E}})}{E\_{\mathcal{P}}} \tag{7}$$

#### *2.6. Statistical Analysis*

The significant role of the independent factors as well as the interactions between these factors was analysed using ANOVA (*p* < 0.05). The contribution of independent factors in the adsorption as a function for the quadratic model was identified by the adjusted coefficient of determination (R<sup>2</sup> adj). The relationship between the independent factors and their effects on the adsorption was represented in a three-dimensional graphical representation, which is RSM.

#### **3. Results and Discussion**

#### *3.1. Characteristics of Durian Adsorbent*

The surface morphology of the durian shell before and after the adsorption process is presented in Figure 2. It was noted that the adsorbent surface having a smooth surface with no pores might be due to the existence of pectin and lignin. These findings are inconsistent with those reported by Lazim et al. [22]. In contrast, the durian shell after the adsorption process have rough surface morphology with dye residues.

The FT-IR spectrum of the durian shell is shown in Figure 3. The FTIR spectra showed a broad peak at 3413.86 cm−<sup>1</sup> which specifies the existence of phenol group of lignin and cellulose. The peaks at 1623.92 are assigned as the C=C aromatic bond or amide areas that are characteristics of proteins and enzymes. The spectrum exhibited at 1437.20 is classified as the C-O stretching vibrations of carboxylic acids (-COOH) and the spectrum at 1054.42 cm−<sup>1</sup> corresponds to the asymmetric stretch of C-N bonds from substituted amines [36]. The small band at 616.98 cm−<sup>1</sup> owing to the C-Cl stretching of alkyl halides [37]. Similarly, Tham et al. [38] observed 3334 cm−<sup>1</sup> which is related to the hydroxyl (O-H) group such as alcohol and phenols. They also observed peaks at 1024 cm−<sup>1</sup> and 1420 cm−<sup>1</sup> respectively representing coupled C–O stretching and amides functional groups of primary amides and aliphatic amides respectively.

**Figure 2.** Surface morphology of durian shell adsorbent using FESEM; (**A**) before (**B**) after.

**Figure 3.** FTIR analysis of durian shell.

#### *3.2. Optimisation of Adsorption of C.I BB16*

The parameter for optimisation of the adsorption of C.I BB16 using durian shell as adsorbent was performed in the current study by Design-Expert software—the maximum decolourisation of C.I BB16 80.7% recorded at pH 4, 30 min, 0.1 g/L of the adsorbent dosage and 10 mg/L of the initial concentration of C.I BB16. In contrast, the maximum reduction of COD was 86.2% at pH 6, after 30 min, with 1 g/L of adsorbent dosage and 10 mg/L of C.I BB16 initial concentrations (Table 1).

*x*<sup>1</sup> (pH), *x*<sup>2</sup> (time), *x*<sup>3</sup> (adsorbent dosage (g/L), *x*<sup>4</sup> (initial concentration of C.I BB16) (mg/L), *y*<sup>1</sup> (decolourisation of C.I BB16) (%), *y*<sup>2</sup> (COD removal) (%).

In comparison Lazim et al. [22] revealed that a durian peel had removed 69.60% of bisphenol with an adsorption capacity of 4.20 mg/g for 24 h. These differences would be related to the types of pollutants treated in their investigations. Lazim et al. [22] studied the bisphenol, while in this study, C.I BB16 dye was treated. Moreover, they treated the durian peel with sulfuric acid before the adsorption bisphenol. Besides, Abd Latib et al. [39] revealed that 99% of methylene blue had been removed using durian shell activated carbon at the conditions 0.60 M H2O2 at 700 ◦C for 30 min. The findings in the present study indicated that the investigated factors play an important role in the adsorption and mineralisation rate of C.I BB16. Decolourisation and mineralisation in the adsorption process are affected by pH, time, durian shell dosage and initial C.I BB16 concentration. In this study, the effect of pH on the adsorption and mineralisation of C.I BB16 was studied at different pH values ranging from 4.0 to 8.0. It is well accepted that the pH of the solution plays a vital part in the adsorption process, which influences the surface charge of the adsorbent and the degree of ionisation and specification of the adsorbate. Generally, it can be deduced that adsorption of C.I BB16 dye by durian shell adsorbent was low at low pH. Therefore, the findings recommended that the adsorption of C.I BB16 on durian shell is favourable in alkaline solutions compared to in neutral and acidic conditions. This is because at the low pH environment, the competition between excess hydroxyl ions, H+, and the cationic groups on the C.I BB16 dye for adsorption sites will make the surface of the adsorbent to exist as a positively charged molecule [40,41].

The adsorption of C.I BB16 investigated the contact time of 30 min to 240 min. It can be deduced that the adsorption and mineralisation rate of C.I BB16 increases as contact time increases. However, a further increase in contact time did not enhance the C.I BB16 decolorisation process and mineralisation rate. The quick absorption rate at the beginning is attributed to the availability of plentiful active binding sites on the adsorbent surface. Nevertheless, the adsorption process develops an attachment-controlled route owing to less accessibility of adsorption sites as contact time increases. Besides, different durian shell dosages were investigated, ranging from 0.1 to 1 g. The results indicated the decolourisation mineralisation rate of C.I BB16 increased when the durian shell dosage increased. This is assigned to a rise in the adsorptive surface area and the accessibility of more active binding sites on the surface of the durian shell. However, the presence of excess durian shell decreases the adsorption capacity of C.I BB16—the decrease in total surface area available for C.I BB16 from aggregation or overlying of adsorption sites [42]. Hence, as the dosage of durian shell adsorbent increases, the amount of C.I BB16 adsorbed onto durian shell adsorbent gets reduced, thus causing a decrease in adsorption capacity. In some cases, a further increase in adsorbent dose did not remarkably alter the adsorption capacity. This is owing to the non-accessibility of active sites on the adsorbent and establishment of equilibrium between the C.I BB16 and the durian shell adsorbent [43]. The rate of adsorption and mineralisation rate is a function of the initial concentration of the adsorbate, which makes it an important factor to be considered for effective adsorption—the percentage removal of C.I BB16 decreased with increasing the initial C.I BB16 concentration (Figures 4 and 5). It seems that all adsorbents have a limited number of active sites and at a certain concentration, the active sites become saturated. The initial concentration may provide the necessary driving force to overcome the mass transfer resistance of C.I BB16 between the aqueous and the solid phase [44]—the increase in the initial C.I BB16 concentration also enhances the interaction between the C.I BB16 anion in the aqueous phase and the biomass surface; resulting in higher uptake of C.I BB16 for the given amount of biomass [29,42].

**Figure 4.** Three-dimensional response surface plot for adsorption of C.I BB16 in aqueous solution;. *y*1: decolourisation; as a response of interaction between independent factors *X*<sup>1</sup> (*pH*), *X*<sup>2</sup> (*Time*), *X*<sup>3</sup> (adsorbent dosage) and *X*<sup>4</sup> (BB16 conc.). (**A**) interaction between X1,X2; (**B**) interaction between X1,X3; (**C**) interaction between X1,X4; (**D**) interaction between X2,X3; (**E**) interaction between X2,X4; (**F**) interaction between X3,X4.

**Figure 5.** Three-dimensional response surface plot for adsorption of C.I BB16 in aqueous solution;. *y*2: COD removal; as a response of interaction between independent factors *X*1, *X*2, *X*<sup>3</sup> and *X*4*X*<sup>1</sup> (*pH*), *X*<sup>2</sup> (*Time*), *X*<sup>3</sup> (adsorbent dosage) and *X*<sup>4</sup> (BB16 conc.). (**A**) interaction between X1,X2; (**B**) interaction between X1,X3; (**C**) interaction between X1,X4; (**D**) interaction between X2,X3; (**E**) interaction between X2,X4; (**F**) interaction between X3,X4.

The regression coefficients for the independent factors are illustrated in Table 2. The interactions between these factors and the linear and quadratic regression coefficients were identified by the least square method. *p*-value with 95% of the confidence level was used as a sign for the importance of each independent factor. The results noted that factors *x*1, *x*2, *x*<sup>3</sup> and *x*<sup>4</sup> have a significant quadratic effect on the decolourisation of C.I BB16 (*p* < 0.0002), as well as for COD removal (%). Moreover, these factors have a synergistic effect on the decolourisation of C.I BB16, final pH and COD removal (%) (*p* < 0.05). The second-order model described the relationship between decolourisation and COD reduction and selected factors as given by final equations in terms of actual factors as following:

$$\begin{array}{rcl} y\_1 &=& 57.76 - 1.09 \mathbf{x}\_1 - 12.41 \mathbf{x}\_2 + 7.73 \mathbf{x}\_3 + 0.20 \mathbf{x}\_4 + 7.55 \mathbf{x}\_1^2 A^2 + 3.90 \mathbf{x}\_2^2 - 7.97 \mathbf{x}\_3^2 \\ &+ 4.20 \mathbf{x}\_4^2 - 0.26 \mathbf{x}\_1 \mathbf{x}\_2 - 1.63 \mathbf{x}\_1 \mathbf{x}\_3 + 1.32 \mathbf{x}\_1 \mathbf{x}\_4 + 4.65 \mathbf{x}\_2 \mathbf{x}\_3 + 0.36 \mathbf{x}\_2 \mathbf{x}\_4 \\ &+ 2.83 \mathbf{x}\_3 \mathbf{x}\_4 \end{array} \tag{8}$$

$$\begin{array}{l} y\_2 = & 35.79 + 2.130\mathbf{x}\_1 - 14.64\mathbf{x}\_2 + 15.55\mathbf{x}\_3 - 6.73\mathbf{x}\_4 - 2.66\mathbf{x}\_1^2 + 16.90\mathbf{x}\_2^2 + 3.31\mathbf{x}\_3^2\\ &+ 1.41\mathbf{x}\_4^2 + 0.87\mathbf{x}\_1\mathbf{x}\_2 - 1.42\mathbf{x}\_1\mathbf{x}\_3 + 5.00\mathbf{x}\_1\mathbf{x}\_4 + 4.65\mathbf{x}\_2\mathbf{x}\_3 + 1.37\mathbf{x}\_2\mathbf{x}\_4\\ &- 0.48\mathbf{x}\_3\mathbf{x}\_4 \end{array} \tag{9}$$

**Table 2.** Regression coefficient and their significance of the quadratic model for the adsorption of C.I BB16 in aqueous solution.

*x*<sup>1</sup> (pH), *x*<sup>2</sup> (time), *x*<sup>3</sup> (adsorbent dosage (g/L), *x*<sup>4</sup> (initial concentration of C.I BB16), *y*<sup>1</sup> (decolourisation of C.I BB16) (%), *y*2, (COD removal) (%).

The summary of the analysis of variance (ANOVA) for the quadratic model is shown in Table 3. It can be deduced that the linear model for decolourisation (%) and COD removal (%) was significant at a confidence level of (*p* < 0.05) with determination coefficients (R2 adj) equal to R2 (adj) are 0.8879 and 0.8603 respectively.

**Table 3.** Analysis of the variance (ANOVA) of the response surface quadratic model for the adsorption of C.I BB16 in aqueous solution.

*y*<sup>1</sup> (decolourisation of C.I BB16) (%), *y*<sup>2</sup> (COD removal) (%); *y*<sup>1</sup> R<sup>2</sup> = 0.9408; R2 (adj) = 0.8579; *y* R<sup>2</sup> = 0.9418; R2 (adj) = 0.8603.

These results ascertain the suitability of the model. The lack of fit for the model was non-significant value (*p* > 0.05), which specifies that the model was predicted for the decolourisation (%) and COD removal (%), while was suitable for prediction within the range of the independent factors explored in the current investigation. Moreover, the small residual values are an emblem for the worthy agreement of the experimental data with the mathematical model. RSM with CCD analysis was executed to examine the interactions between the variable factors within the tested range when one factor is an actual factor at the centre point. For the decolourisation (%) and COD removal (%), a significant (*p* < 0.05) interaction between *x*<sup>1</sup> (4–8), *x*<sup>2</sup> (30–240 min), *x*<sup>3</sup> (0.1–1 g) and *x*4(10–20 ppm) as actual factors is shown in Figures 4 and 5. The correlation plot between model-predicted value vs. actual value and residuals vs. run number; A; decolourisation (*y*1), B; COD removal (*y*2) is depicted in Figure 6. The results revealed a significant correlation between the independent factors and the removal percentage of the dependent variables.

**Figure 6.** Correlation plot between model-predicted value vs. actual value and residuals vs. run number; (**A**) decolourisation (*y*1), (**B**) COD removal (*y*2).

Validation of the Optimal Parameters

The possible direction for maximising the adsorption of C.I BB16 is by using point optimisation technique with expert design software. The optimal operation parameters of the adsorption of C.I BB16 were recorded at pH 8 (*x*1), 30 min (*x*2), 1.0 g of adsorbent (*x*3) and 15 ppm of C.I BB16 (*x*4). At the optimal conditions of the adsorption of C.I BB16, the independent factors interacted significantly (*p* < 0.05) (Table 4).

**Table 4.** The best operating parameters for the adsorption of C.I BB16 in aqueous solution.

*x*<sup>1</sup> (pH), *x*<sup>2</sup> (time), *x*<sup>3</sup> (adsorbent dosage (g/L), *x*<sup>4</sup> (initial concentration of C.I BB16), *y*<sup>1</sup> (Decolourisation of C.I BB16) (%), *y*<sup>2</sup> (COD removal) (%).

This specifies that interaction effects occur as a consequence of one factor influenced by another factor. The pH, time, durian shell adsorbent dosage and initial concentration of C.I BB16 act together very well to attain the maximum decolourisation and COD removal rate. The percentage error among the actual and predicted values was recorded less than 5% (confidence level 95%) which indicates that no substantial difference was observed and therefore, affirms that the RSM method was suitable to optimise the operational settings of adsorption of C.I BB16. In comparison with the previous studies by Mohammed et al. [19], it was noted that durian shell sorbents achieved 95.91% and 97.81% decolourisation of C. I Methylene Blue and C. I Brilliant Green, respectively. The differences in the results occurred as these factors were investigated individually, whereas, in their study, the interactions between these independent factors were not considered. In the current investigations, the optimisation of the adsorption of C.I BB16 achieved more than 77.61%, which designate the role of RSM design in attaining the maximum and optimum conditions for the adsorption of C.I BB16. The high effectiveness of adsorption of C.I BB16 in reaction to the examined independent variables would be described based on the occurrence of a significant synergic relationship between the parameters. These outcomes are constant with the authors in the literature, who directed that the investigated factors play a significant part in the adsorption of C.I BB16, Mondal et al. [29]. The optimisation process using RSM regulates the best operation parameters to reach high adsorption efficiency.

#### *3.3. Mechanism of Adsorption*

To describe the nature of adsorption, the mechanism of adsorption was studied by interpretation of the intraparticle diffusion (Figure 7).

The key factors that control the rate of adsorption can be affected by the adsorption at the internal and external surface which depends on the binding process, transfer of mass through the external boundary and the adsorbate diffusion to the adsorption site across liquid pores or solid surface. Therefore, the decolourisation (77.61%) and COD removal (80.60%) at pH 8, 30 min, 1 g/L of durian shell dosage and 15 ppm of C.I BB16 were achieved based on independent factors explored in the current investigation. The competition between hydroxyl ions, H+ and the cationic groups on the C.I BB16 dye for adsorption sites will change the surface of adsorbent charged molecule positively or negatively depends on pH value, which will improve the positively charged dye cations through electrostatic forces of attraction. Durian is the member of the Bombacaceae family which is largely found in Malaysia and Indonesia of Southeast Asia. Durian peels are the waste materials of durian and consist of a lot of polysaccharides, cellulose and lignin. Those ingredients have different active functional groups that relatively affect the adsorption of analyses (dye, metal), via interactions occurred between solid surfaces and adsorbate (Ahmad and Danish, 2018); structure of the adsorbate and the functional groups of the adsorbent surfaces. Durian adsorption of C.I BB16 molecules results from hydrogen bonding between OH and NH2 groups of dye molecules and OH groups of durian peels and electrostatic interaction between the negatively charged durian surface and cationic dye. Referred to FTIR spectra, the results showed that carboxylic acids (-COOH) and the hydroxyl groups were available on the sites because of interactions of the surface of the peel. These functional groups could control the uptake of positively charged C.I BB16 molecule. Whereas, electrostatic attraction, h-bonding could occur between the negatively charged surface and positively charged +NH-CH3-C.I BB16 molecules [45,46].

**Figure 7.** Removal mechanism of C.I BB16 using durian shell adsorbent.

*3.4. Techno-Economic Analysis for Preparation and Application of Durian Shell for Hair Dye Greywater Treatment*

The durian shell adsorbent preparation and application diagram was designed using SuperPro Designer, as presented in Figure 8b. The suggested adsorbent preparation and application unit is proposed to have consisted of four stages. In stage A, the raw materials of durian shell adsorbent are collected, washed and granulated as well as adsorbent dose preparation is done. In stage B, the hair dye greywater is subjected to primary treatment to remove the large suspended solids, which might affect the adsorption process negatively. The main treatment process of the hair dye greywater by the adsorption using durian shell adsorbent is supposed to be performed in Stage C. Finally, the stage proposes the post-treatment process for the treated hair dye wastewater which is safe for disposal while the durian shell might be used as the substrate for fuel or biohydrogen production because of the high contents of the carbohydrates.

The FCE represents all equipment purchase cost, process piping, equipment installation, electrical systems, instrumentation and controls, buildings, construction and yard improvements as well as the *WCC* which might be up to 6.5% of FCE as described by Han et al. [46]. Therefore, the FCE for designing a preparation and coagulation unit for wastewater treatment with 1000 m3/day capacity could be up to USD 32,000.00 (Table 5).

**Figure 8.** (**a)** The flowchart for preparation and application of durian shell adsorbent for C.I BB16 wastewater treatment. (**b**) Internal rate of return (IRR) of durian shell adsorbent production.

The additional costs for the treatment process might reach 50% of the equipment cost, which was equivalent to USD 16,000.00 [47]. The TCI of the plant might go up to USD 48,000.00 WCC, which was estimated to be 6.5% of the FCE to the tune of about USD 3120.00 Therefore, according to Equation (1), the TCI of USD 51,120.00 could be calculated. The techno-economic analysis is one of the best methods for determining the effectiveness of any alternative method for treating environmental pollution because this process provides more details for the total cost of the treatment process for specific pollution in comparison to the current technologies [47]. It has also been used to estimate the H2 production by fungi from food wastes and provided very accurate information on the total cost [44,48].

**Table 5.** Fixed capital estimate (FCE) for preparation and application of durian shell adsorbent for hair dye wastewater treatment.

\* Refer the item code to Figure 8a.

#### 3.4.1. Annual Operation Cost

The *CRM* is the cost of durian shell. The durian shell as well as the chemicals required for preparation of the durian shell to be used as a coagulant was supplied by a local supplier. The *CU* include the cost of electricity and water which are required for the preparation of durian shell adsorbent and the wastewater treatment and was estimated based on the price for each unit in the local currency. The *CWG* represents the final solid waste generated from the coagulation treatment unit of the wastewater which is considered as a biomass yield with high concentrations of carbohydrates and could be used as bio-sorbents for wastewater treatment as well as for fuel production. Moreover, the generated wastewater from the adsorption unit could be safe for disposal.

Table 6 shows the annual operation cost for the preparation of durian shell adsorbent. The raw material cost (*CRM*) was USD zero/year (since the durian shells are agriculture wastes and no chemical are required for the preparation). In the current study, the best dose for achieving high COS, turbidity and TSS removal are 1 g/L. The suggested design here manages with 1000 m3 per day for 240 days per year. The experiments confirmed that 1 kg of raw durian shell generates 500 g pure coagulant. In order to deal with 240,000 m3 per year, the process should deal with 480 tons of raw durian shell to generate 240 tons (1000 kg/day) of pure coagulant to treat the quantities required per year. The total cost of utilities is estimated to be USD 5000.00/year as estimated by Han et al. [48] in China, which has a similar cost for water and electricity as Malaysia. The costs of maintenance, insurance and labour as well as the operating labour cost, were based on the estimation that each process operation using 100 kg of raw durian shell to produce 50 kg of pure adsorbent daily required only one operator with an average salary of 6000.00 USD/year per operator (max 3 h/day). The maintenance and insurance were estimated at 2% and 1% of the FCC, respectively, according to Vlysidis et al. [47] and Ljunggren and Zacchi [49]. It can be concluded that the AOC is USD 11,960.00/year.

**Table 6.** Annual operation cost of the preparation and application of durian shell adsorbent for hair dye wastewater treatment.

#### 3.4.2. Annual Profitability and Revenue

The revenue for the preparation and application of durian shell adsorbent is generated from selling the adsorbent for the treatment of different wastewater. The current price of chemical adsorbent might reach USD 261.81.00 per tons of Honeydew peels-activated carbon (HDP-AC) and the commercial market of the activated carbon price is USD 1000.00/tons [50]. The machines in Table 5, have 100 kg of the capacity for each operation which takes place within a day. Therefore, the total capacity per year is to produce 240 tons of pure durian shell adsorbent. The total net of the raw durian shell in the current work was 500 g per 1000 g of the raw materials. Therefore, the revenue for each production process is 1000 kg of the pure durian shell adsorbent which gives USD 62834.00 of the revenue. The annual revenue of USD 62,834.00 could be estimated. Considering the 15% local tax in Malaysia, this is estimated to be USD 9425.00 and the annual operation cost is USD 11,960.00, the annual profitability (after the tax) of the durian shell adsorbent is USD 41,499.24 (Table 7a). Based on Equation (9), the specific cost of pure durian shell adsorbent production is estimated to be USD 172.71 per ton which is lower than the market price of honeydew peels-activated carbon (HDP-AC) (USD 261.81) and the commercial market of activated carbon price which is USD 1000.00/tons.

**Table 7.** (**a**) Annual profitability of the preparation and application of durian shell adsorbent for hair dye wastewater treatment. (**b**) Cost factors associated with the adsorption of C.I BB16 from the hair dye wastewater by durian shell adsorbent for hair dye wastewater in comparison to the physical and chemical treatment process.

In this study, the pure durian shell adsorbent production unit was economically evaluated according to the internal rate of return (IRR), net present value (NPV) as well as the payback period (PBP) of 10 years of the lifetime. The IRR (the efficiency of the investment) was 45.62% in the current study (Figure 8b). Based on the analysis given above, the durian shell exhibited high efficiency for adsorption of hair dye wastewater. Therefore, the techno-economic analysis as described in the current work is economically feasible by 45.62% of the IRR. It might provide alternative methods for chemical and physical coagulation for wastewater treatment (Table 7b).

#### **4. Conclusions**

It can be concluded that an adsorbent was successfully produced via the durian shell and successfully showed high efficiency for the adsorption of C.I BB16—the optimal operating factors for adsorption C.I BB16 recorded at pH 8, time (30 min), durian shell dosage (1 g/L), and 15 mg /L of C.I BB16 concentrations where 77.61 vs. 74.26 (%) of C.I BB16 removal and 80.60 vs. 78.72 (%) of COD removal was the observed and predicted results recorded with an R<sup>2</sup> coefficient of 0.94, respectively. The process is very economically feasible with 45.62% of the internal rate of return (IRR). These findings indicate that the durian shell adsorbent has great potential to be applied as a low cost and eco-friendly adsorbent in the adsorption process of dyes in aqueous solution. The utilisation of durian shells as an adsorbent also helped in reducing solid waste in agriculture.

**Author Contributions:** Conceptualization, A.A.-G.; methodology, Y.G. and M.M.A.; software, M.A.-S.; validation, A.A.-G., M.A.M. and R.M.S.R.M.; formal analysis, A.A.-G.; investigation, S.A.; resources, A.A.-G.; data curation, Y.G.; writing—E.N., A.A.-G.; writing—review and editing, M.A.M.; visualization, A.A.-G.; supervision, R.M.S.R.M.; project administration, M.A.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was funded by the technical & financial support from UNITEN RMC Internal Research Grant: RJO 10517919/iRMC/Publication and Ministry of Education Malaysia (KPM) through the Fundamental Research Grant Scheme (FRGS) with reference code: FRGS/1/2019/WAB05/UTHM/02/5 (Modification of Bead Adsorbents with Ceramic Sanitary Ware Waste (CSWW) and Chitosan for Laundry Greywater (LGW) Safe Disposal) as financial support for this research project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Diversity of Carbon Storage Economics in Fertile Boreal Spruce (***Picea Abies***) Estates**

**Petri P. Kärenlampi**

Lehtoi Research, 81235 Lehtoi, Finland<EMAIL_ADDRESS>
**Abstract:** A "normal forest", an idealized estate with a uniform distribution of stand ages, can be used in the study of sustainable management practices. As the normal forest contains a variety of stand ages, the characteristics of the stands can be represented in terms of a "normal stand", with properties known as a function of age. This paper takes seven never-thinned stands as seven "normal stands", which describe seven estates of normal forest. The intention is to study the robustness of carbon storage microeconomics to varying estate characteristics. It was found that the economically optimal rotation ages vary. The state sums of volume and capitalization, corresponding to any optimal rotation, also vary significantly. Growth rates vary more than the optimal expected stand volumes. Consequently, any excess volume related to carbon storage adds on to an almost unified basic volume. For all seven normal estates, the most economical way of increasing carbon storage is to increase the size of trees retained in thinning from above.

**Keywords:** capitalization; capital return rate deficiency; expected value; carbon storage; estate fertility; carbon rent

#### **1. Introduction**

Boreal forests constitute a significant carbon sink. A particular benefit of boreal regions is carbon storage in the soil. The amount of soil carbon may exceed the carbon storage in living biomass [1–4]. However, living biomass produces the litter that results in soil carbon accumulation, and consequently the rate of carbon storage depends on the rate of biomass production on the site. The biomass production rate is related to the amount of living biomass [2,5–7].

In the occurrence of clearcutting and soil preparation, the net release of carbon from the soil to the atmosphere begins [8–11]. However, recent unpublished datasets indicate that maintenance of canopy cover may or may not be essential in carbon sequestration.

Human activity related to the sequestration and release of carbon into the atmosphere forms a complex system of industries and disciplines [12–16]. Instead of trying to propose any interdisciplinary optimum, we here discuss microeconomics of carbon sequestration in forestry, focusing on the estate level.

In principle, it would be straightforward to apply a carbon trade policy in forestry, possibly administered by public agents [17]. An unbiased carbon trade, however, would incorporate a high initial expense [18]. However, it has been recently shown that a carbon rental system is equivalent to the carbon trade system, without much of an initial expense [18].

From the viewpoint of generic instructions or policy actions, it might be beneficial to reduce the variety of initial estate states by adopting some kind of unifying boundary conditions. A tempting candidate is the normal forest principle [19]. This principle simply refers to the postulation that stand ages are evenly distributed, and stand characteristics are uniquely determined by stand age. Such a postulation, even if often departing from reality, simplifies many treatments, producing idealized systems that are stationary in time. Application of the normal forest principle allows for the determination of an optimal

**Citation:** Kärenlampi, P.P. Diversity of Carbon Storage Economics in Fertile Boreal Spruce (*Picea Abies*) Estates. *Sustainability* **2021**, *13*, 560. https://doi.org/10.3390/su13020560

Received: 15 December 2020 Accepted: 6 January 2021 Published: 8 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

rotation age, as well as expected values of estate characteristics corresponding to the financial optimum.

Microeconomically optimal rotation age and the expected value of stand volume are not necessarily optimal from the viewpoint of national economics or the viewpoint of carbon sequestration. It has been recently shown that microeconomics often favors solutions with relatively low capitalization [20–23]. Low capitalization reduces volumetric growth rate, as well as litter accumulation rate [21,23,24].

A rather recent paper clarified the microeconomics of carbon sequestration within a boreal spruce forest estate [23]. It found that the best capital return rate is achieved by multiple repeated thinning from above to the transition diameter where sawlogs are gained instead of pulpwood only [21–23]. Increased carbon storage was found to be most economically achieved by increasing the harvesting limit diameter [23]. However, carbon rent derivable from present carbon market prices appeared not high enough to compensate for the induced capital return deficiency [23].

A significant revision to the recent approach [23] has been presented still more recently [25]. Literature values of the yield of sawlogs and pulpwood [26,27] were replaced by empirical observations [25]. The resulting value increment of a typical spruce tree as a function of breast-height diameter is a much smoother function than the literature values would indicate [25]. Some small trees have yielded sawlogs. On the other hand, the expected sawlog content of larger trees is smaller than the literature values would indicate, as operators are trained not to produce sawlogs that may be rejected at the mill [25].

The consequences of smoother value increments are significant [25]. Due to the price premium of sawlogs from clearcutting, all feasible silvicultural treatments on fertile spruce estates terminate in clearcutting instead of multiple repeated thinnings from above [21–23,25]. Financially optimal expected values of timber stocks per hectare are greater when the empirical log yield function is used, in comparison to the literature values of sawlog yield [21,23,25].

The recent treatments have been based on field observations, a growth model, and financial theory. The growth model, as well as the financial treatment, apply probability theory. However, the outcome of the growth model is deterministic, since it is composed in terms of expected values only [21,22,28]. The financial treatment, even if probabilistic as such, applies the normal forest principle, thus applying the postulate of constant stand age distribution within an estate [21–23,25,29,30]. Consequently, any density function of stand age becomes determined by one parameter (rotation age), and the financial treatment becomes deterministic.

The field observations naturally contain statistical scatter. However, recent papers extended the normal forest principle in a somewhat novel manner [23,25]. Instead of discussing the variety of stand properties in detail, a representative "normal stand" was composed as a combination of two stands of average fertility and younger-than-average age. Application of this normal stand procedure has effectively neglected statistical scatter in field observations and rendered the entire treatment deterministic [23,25].

This paper intends to clarify the effects of statistical scatter in field observations. It is of interest whether the results gained when such a scatter is neglected are robust in the presence of the scatter. However, we do not abandon the normal stand - principle. Instead, we applied that principle separately for any of the seven experimental stands where field data has been collected. Each of the seven unthinned stands is used as a normal stand characterizing an estate. The seven normal forest estates is then discussed in order to clarify the robustness of recent results to statistical scatter in field observations. It is worth noting that discussion of any of the seven constituted estates is deterministic—it is of interest whether some of them significantly differ from results gained by neglecting the scatter in field observations.

First, experimental materials are presented, as well as related computational methods. Then, microeconomic methods, as well as carbon rent formulae are introduced. The methods largely unify with a very recent paper [25], apart from the consideration of the scatter

in the characteristics of the normal stands. Third, results are reported for four different kinds of treatment schedules. Finally, the outcome is discussed, and conclusions regarding the robustness of carbon sequestration systems to estate characteristics are given.

#### **2. Materials and Methods**

#### *2.1. Empirical Observations*

Eleven circular plots of 314 square meters in area were taken in November 2018 from typical spots of 11 spruce-dominated forest stands in Vihtari, Eastern Finland. The area of the Vihtari estate is 30 ha, with an elevation from 115 to 125 m above sea level. Seven of the stands had experienced only young stand cleaning, whereas four of the stands were previously thinned commercially. Breast-height diameters were recorded as well as tree species, and a quality class was visually determined for any measured tree, reflecting its suitability for adding value in the future.

On measured plots on sites without any previous commercial thinning, the basal area of tree trunks at breast-height varied from 32 to 48 m2/ha, and the stem count from 1655 to 2451 per hectare. On measured plots on stands previously thinned commercially, the basal area varied from 29 to 49 m2/ha, and the stem count from 891 to 955 per hectare. Further details of the experimental stands are reported in an earlier paper [21].

The Vihtari estate was characterized by measurements from the 11 sample plots. However, such a sample did not necessarily represent the entire estate accurately. More importantly, the sampling did not conform to the normal forest principle, with assumptions of constant age distribution and stand characteristics uniquely determined by stand age. Here, we utilize the normal forest principle in terms of establishing a "normal stand" based on the observations and then approximating the development of this "normal stand" as a function of stand age. Introducing an even stand age distribution produces a "normal estate" on the basis of the normal stand, following the normal forest principle [20].

In recent studies [23,25], the normal stand has been established based on two neverthinned sample plots of medium site fertility and a younger range of stand age. The intention was to establish one normal stand representing the typical characteristics of the 11 sample plots measured at the Vihtari estate. The selection of such a representative stand simplified the treatments and rendered the entire analysis deterministic. In other words, it neglected the variation between the 11 sample plots.

The intention of this paper being to investigate the variability among the sample plots, any of the seven never-thinned sample plots is taken as a normal stand. Commerciallythinned stands are neglected since the thinning revenue i not known. Any of the normal stands iss taken to represent a normal estate, representing a "normal forest" [20]. The robustness of recent results to the scatter among sample plots is then investigated in terms of the seven normal estates.

Instead of using literature values [21–23,26] as an estimate of the yield of sawlogs and pulpwood from any tree of particular breast-height diameter, a recently collected dataset of 6123 spruce trees is used [25]. The data was collected by four different single-grip harvesters operated by six individuals at four harvesting sites. Log cutting instructions provided by three different sawmilling companies were applied. One thinning site was located at Vihtari, and two thinning sites and one clearcutting site at Ilomantsi, all in eastern Finland. The area of the sampling territory was 100 km (west to east) by 30 km (south to north). Elevation varied from 115 to 205 m above sea level.

The proportion of sawlogs of the total harvester-measured volume, as a function of breast-height diameter, is shown in Figure 1 [25]. Figure 1 also shows the literature values used in [21–23,26]. It is found that the empirical sawlog yield function is smoother than the function taken from the literature. Some sawlogs were gained from trees of a diameter of 175 mm, unlike in the literature data. On the other hand, the sawlog proportion of trees thicker than 200 mm is smaller in the empirical function. The latter probably is due to harvester operators being trained not to produce sawlogs that would be rejected at the mill.

**Figure 1.** Literature values as well as empirical values of the sawlog proportion within the commercial section of spruce stems of different sizes.

#### *2.2. Growth Model*

For prognostication of further development of the normal stand, some kind of growth model is needed. The growth model of Bollandsås et al. [21,22,28] was adopted, discussing not only growth but also mortality and recruitment. The original growth model [21,22,28] discussed 50 mm breast-height diameter classes within a temporal resolution of five years. Any diameter class was represented by its central tree, and the process of growth was described in terms of the probability of any tree to transfer to the next diameter class [21,22,28]. The underlying idea of the description of growth is that any tree either remains in the same diameter class or transfers to the next diameter class within the five-year time interval. This underlying idea naturally greatly simplifies computation.

It's found from Figure 1 that according to the literature values of sawlog yield, there is a sharp transition between diameter classes centered at 175 and 225 mm. Such a sharp transition would induce a huge value increment. With the empirical yield function being smoother, there is a need for a greater resolution in the tree size description. It is not too complicated to modify the growth model from the size resolution of 50 mm to 25 mm [25]. However, to retain the underlying principle of the growth model, this requires a corresponding change in the temporal resolution, from 5-year time steps to 30 months.

The simultaneous change in the time and size resolutions retains the probability of any tree to transfer to the next diameter class. The same is not true in the case of recruitment and mortality. The latter two quantities are affected by time resolution only. Correspondingly, the recruitment and mortality values become scaled along with the time step by a factor of 1/2.

For any 25-mm diameter class of trees, the volumetric amount of two assortments, pulpwood and sawlogs are taken from the empirical observations as the expected values for any diameter class.

Description of stand development until the time of observation in 2018 required another kind of approach. An exponential growth function was fitted between an approximated worthless initial volume of 15 m3/ha and the 2018 commercial volume estimate for any of the seven sample plots.

Possibly, the simplest way to approximate financial history is to determine an internal rate of operative return for the period from stand establishment to 2018. In other words, we require

$$V(\tau\_3)e^{-s\tau\_3} - R(\tau\_1)e^{-s\tau\_1} - C(\tau\_2)e^{-s\tau\_2} = 0\tag{1}$$

where *R* is regeneration expense at regeneration time *τ*1, *C* is young stand cleaning expense at cleaning time *τ*2, and *V* is stumpage value of trees at observation time *τ*3. In this study, the observation time *τ*<sup>3</sup> corresponds to November 2018, regeneration time is clarified according to known stand age, and young stand cleaning is assumed to have occurred ten years after regeneration. It is assumed that prices and expenses do not evolve in real terms, and thus presently valid expenses can be used in Equation (1). The regeneration expense is taken as 1250 Euro/ha, and young stand cleaning as 625 Euro/ha.

It is worth noting that the operative internal rate of return *s* in Equation (1) does not correspond to the capital return rate in the entire activity, as the latter depends on nonoperative capitalization such as bare land value.

#### *2.3. Financial Considerations*

To determine a momentary capital return rate, we need to discuss the amount of financial resources occupied [21,23,24,31]. This is done in terms of a financial potential function, defined in terms of capitalization per unit area K. The momentary capital return rate becomes

$$r(t) = \frac{d\kappa}{K(t)dt} \tag{2}$$

where *κ* in the numerator considers value growth, operative expenses, interests, and amortizations but neglects investments and withdrawals. In other words, it is the change of capitalization on an economic profit/loss basis. *K* in the denominator gives capitalization on a balance sheet basis, being directly affected by any investment or withdrawal. It is worth noting that timber sales do not enter the numerator of Equation (2); selling trees at market price levels does not change the amount of wealth, it only converts wealth from trees into the form of cash. However, harvesting naturally changes capitalization appearing in the denominator of Equation (2). Harvesting also likely changes the change rate of capitalization occurring after the harvest.

Equation (2) gives a momentary capital return rate, not necessarily sufficient for management considerations. By definition, the expected value of capitalization per unit area is

$$\langle K \rangle = \int\_{-\infty}^{\infty} p(K) KdK \tag{3}$$

where *p*(*K*) is the probability density function of capitalization *K*. By change of variables we get

$$
\langle K \rangle = \int\_0^\tau p(K) K \frac{dK}{da} da = \int\_0^\tau p(a) K(a) da \tag{4}
$$

where *a* is stand age (or time elapsed since latest regeneration harvesting), and *τ* is rotation age. The expected value of the change rate of capitalization is

$$
\left< \frac{d\kappa}{dt} \right> = \int\_0^\tau p(a) \frac{d\kappa(a)}{dt} da \tag{5}
$$

Correspondingly, the expected momentary rate of relative capital return is

$$\langle r(t)\rangle = \frac{\left\langle \frac{d\mathbf{x}}{dt} \right\rangle}{\langle K\rangle} = \frac{\int\_0^\tau p(a) \frac{d\mathbf{x}(a,t)}{dt} da}{\int\_0^\tau p(a)K(a,t) da} = \frac{\int\_0^\tau p(a)K(a,t)r(a,t) da}{\int\_0^\tau p(a)K(a,t) da} \tag{6}$$

We find from Equation (6) that the expected value of capital return rate within an estate generally evolves in time as the probability density of stand ages evolves. However, Equation (6) can be simplified to be independent of time by adopting the normal forest

principle, where stand age probability density is constant [20]. Besides, the constancy of the expected value of capital return rate in time requires that prices and expenses do not evolve in real terms. Then, the expected value of the capital return rate becomes

$$\langle r \rangle = \frac{\int\_0^{\frac{\pi}{4}} \frac{d\kappa(a)}{dt} da}{\int\_0^{\frac{\pi}{4}} K(a) da} = \frac{\int\_0^{\frac{\pi}{4}} K(a) r(a) da}{\int\_0^{\frac{\pi}{4}} K(a) da} \tag{7}$$

It has been recently shown that Equation (7) corresponds to the ratio of the partition functions of the change rate of capitalization and capitalization itself [32]. It also has been recently shown that the maximization of the net present value of future revenues may result in financially devastating consequences [33]. Momentary capital return rate as given in Equation (2) was introduced in 1860 [31]; an expected value was mentioned in 1967 [29,30]; however, applications have been introduced only recently [21–23,25,32,33].

In the context of the improvement harvesting, 20% of the stem count of good-quality trees is removed in all diameter classes due to the establishment of striproads.

The stumpage value is determined in terms of roadside price, deduced by harvesting expense. We here use the roadside prices recently applied by Parkatti et al. [34], 34.04 Euro/m3 for spruce pulpwood, and 58.44 Euro/m<sup>3</sup> for sawlogs. We further use the same harvest-expense function as Parkatti et al. [34], stated to be based on a productivity study of Nurminen et al. [35]. In addition to the expense function, we include a fixed harvesting entry expense per hectare. The justification of the latter is that some sites require at least partial pre-harvest cleaning. The entry expense is approximated as 200 Euro/ha. Again, it is assumed that prices and expenses do not evolve in time. In other words, the capital return rate is discussed in real terms.

The original version of the growth model operates in five-year time steps [21–23,28], as discussed above. Consequently, an eventual harvesting entry may take place every five years. According to recent investigations, it often is favorable to harvest every five years [21,23]. However, the fixed harvesting entry cost of 200 Euro/ha restricts low-yield harvesting entries. Numerical investigations indicated that it is not reasonable to harvest if the yield would be less than 20–30 m2/ha. Such an entry limit is in concert with practices applied in the area.

The feasibility of harvesting is investigated in five-year intervals, always considering two options: thinning or clearcut. Clearcutting expenses are lower than thinning harvesting costs, according to Parkatti et al. [34], stated to be based on a productivity study of Nurminen et al. [35]. Moreover, a 15% clearcutting premium for the roadside price of sawlogs is applied, following local tradition. The premium, as well as the harvesting cost reduction, is applied within the last 30-month period before eventual clearcutting.

Thinning procedures are iteratively designed to maximize the expected value of capital return rate, up to the rotation age providing the maximum expected value. After the maximum expected value, any thinning procedure is designed to maximize the capital return rate within the next five years.

Four alternative thinning strategies are investigated for any of the seven normal estates. Two of the strategies are based on thinning from above. The first does not apply any restriction to cutting limit diameters, but the thinning procedure is designed in order to maximize the expected value of capital return rate according to Equation (7). The second does apply a minimal 238 mm diameter limit for thinning from above as a restriction Ref. [25]. The remaining two thinning strategies refer to proportional thinning. In other words, the probability of tree removal is the same in any diameter class. The first of the two proportional thinning strategies is taken as heavy proportional thinning, where 50% of acceptable-quality trees are removed. Gentle proportional thinning removes 30% of acceptable-quality trees in any diameter class.

Application of Equation (7) results in an expected value of capital return rate for any treatment schedule investigated. Treatment schedules not corresponding to the most economical one induce some amount of capital return rate deficiency, in terms of percentage per annum. Such deficiency is related either to the amount of timber stock deviating from the optimal timber stock, or the rotation age deviating from the optimal rotation age, or both. In the former case, the deficiency often is due to an excess standing volume, which allows for expressing the deficiency per excess volume unit.

#### *2.4. Carbon Rent Considerations*

The last issue in this section of methods regards carbon trade and carbon renting. It is worth noting that carbon prices or rents do not enter the analysis described above in any way. Carbon prices and rents are discussed to enable a comparison of the capital return deficiency per excess volume with any hypothetical carbon rent.

It has been recently shown that policies based on carbon rent are equivalent to policies based on carbon sequestration subsidies and taxies [18]. Unbiased carbon sequestration trade would require a huge initial investment; correspondingly, mostly carbon rent procedures are practically feasible [18]. We here present a brief derivation of the equivalency of the two principles of subsidies, although adopting boundary conditions possibly less restrictive than those of Lintunen et al. [18].

Let us establish a carbon sequestration subsidy system at a particular time *τ*1. Within a time range up to time *τ*2, the total carbon trade compensation is

$$p\_{\tau\_1} \mathbb{C}\_{\tau\_1} + \int\_{\tau\_1}^{\tau\_2} \frac{d(p\mathbb{C})}{dt} dt = p\_{\tau\_1} \mathbb{C}\_{\tau\_1} + \int\_{\tau\_1}^{\tau\_2} \left[ p \frac{d\mathbb{C}}{dt} + \mathbb{C} \frac{dp}{dt} \right] dt \tag{8}$$

where *pt* is carbon price at time *t*, and *Ct* is carbon inventory at time *t*. On the other hand, the revenue from carbon rentals is

$$\int\_{1}^{2} uCdt\tag{9}$$

where u is the rent rate per carbon unit. Now, to establish equivalency between the carbon storage trade and rent, Equations (8) and (9) must become equal. Equality naturally should apply in any possible circumstance. One of the circumstances is that the time change rate of prices, as well as inventories, is zero. In such a case, the latter term of Equation (8) vanishes. Consequently, the long-term flow of carbon rents should equal a one-time initial storage purchase payment. If the duration of the rent payments extends towards infinity, the only possibility is that the present value of rent payments forms a contracting series. One possibility of such a contracting series is

$$\mu\_{\mathbb{T}\_1} \mathbb{C}\_{\mathbb{T}\_1} \int\_0^\infty e^{-qt} dt = p\_{\mathbb{T}\_1} \mathbb{C}\_{\mathbb{T}\_1} \tag{10}$$

where *q* is a discount rate. The corresponding solution for the carbon rent rate is

$$u\_t = qp\_t\tag{11}$$

One can readily show that Equation (11) applies not only to a steady-state of Equations (8) and (9) but also to any incremental carbon price and inventory.

#### **3. Results**

The expected value of capital return rate as a function of rotation age, integrated over the stand lifespan according to Equation (7) is shown in Figure 2. In Figure 2, no restriction is applied in the maximization or the capital return rate. Any feasible treatment cycle terminates in clearcutting, and thus utilizes the premium in clearcutting price, as well as the harvesting expense reduction. In most cases, thinnings from the above are made to 188 mm cutting diameter limit, even if also 213 mm limits appear. The number of thinnings up to the financially optimal rotation age varies from one to three. After the optimal rotation age, one thinning is triggered in five of the seven cases, in the absence of prior clearcutting.

**Figure 2.** Capital return rate according to Equation (7) for seven different normal estates, without any restriction in the design of the applied thinning procedure. The legend identifies the experimental plot identifying any normal estate.

It is found from Figure 2 that the normal estates, with properties determined from the corresponding normal stands, differ in terms of the expected value of capital return rate, but they differed more in the optimal rotation age. The amplitudes of variation are 9.8% and 31%, respectively. The optimal rotation ages span from 45 to 85 years. This is due to fertility, as well as the number and size of trees observed. A third significant contributor is the age of any normal stand in the occurrence of data sampling. The oldest stand was of age 45 years. The capital return rate probably would be higher and the corresponding rotation age lower if the stand had been thinned earlier.

The expected value of capital return rate as a function of rotation age is shown in Figure 3, with the restriction of a minimum cutting limit diameter of 238 mm in thinnings from above. Again, any feasible treatment cycle terminates in clearcutting, and thus utilizes the premium in clearcutting price, as well as the harvesting expense reduction. The restriction induces gentler thinnings and greater capitalizations. The restriction naturally reduces achievable capital return rates. The gentler thinnings reduce the optimal rotation age in five of the seven cases, even if it increases the number of thinnings. The number of thinnings up to the financially optimal rotation age varies from one to six. The optimal rotation ages now span from 40 to 85 years, the latter corresponding to the case of delayed first thinning. The amplitude of variation in capital return rate is 10.2%, and 36% in optimal rotation age.

**Figure 3.** Capital return rate according to Equation (7) for seven different normal estates, with a minimum 238 mm cutting diameter limit in thinnings from above. The legend identifies the experimental plot identifying any normal estate.

The expected value of capital return rate as a function of rotation age is shown in Figure 4, with heavy proportional thinning. Again, any feasible treatment cycle terminates in clearcutting and thus utilizes the premium in clearcutting price, as well as the harvesting expense reduction. The heavy proportional thinning renders lower capital return rates in comparison to Figures 2 and 3, as well as shorter rotation times. The optimal rotation ages now vary from 35 to 60 years. The amplitude of variation in capital return rate is 12.9%, and 26.3% in optimal rotation age. In all cases, there is only one thinning. However, in six of the seven cases the financial result is inferior in comparison to no thinning at all.

**Figure 4.** Capital return rate according to Equation (7) for seven different normal estates, with heavy proportional thinning. The legend identifies the experimental plot identifying any normal estate.

The expected value of capital return rate as a function of rotation age is shown in Figure 5, with gentle proportional thinning. Any feasible treatment cycle terminates in clearcutting, and thus utilizes the premium in clearcutting price, as well as the harvesting expense reduction. The gentle proportional thinning renders lower capital return rates in comparison to Figures 2 and 3, as well as shorter rotation times. In comparison to the heavy proportional thinning (Figure 4), capital return rates are higher but optimal rotations shorter. In all cases, there is only one thinning before the rotation age corresponding to the maximum expected capital return rate. Somewhat surprisingly, in five of the seven cases, the financial result is inferior in comparison to no thinning at all. It is, however, worth noting that gaining the clearcutting price premium may be uncertain in the absence of any thinning.

**Figure 5.** Capital return rate according to Equation (7) for seven different normal estates, with gentle proportional thinning. The legend identifies the experimental plot identifying any normal estate.

The capital return rate deficiency in terms of Euros per excess volume is shown in Figure 6, as a function of excess volume, linking economic and environmental sustainability. No restriction is applied in the maximization or the capital return rate. Only positive values are shown in Figure 6. The rotation age providing the maximum capital return rate corresponds to the origin in any series of observations. All the excess volumes in Figure 6 are due to extended rotation age. It is found that even excess volumes of 10 to 25 cubic meters per hectare induced a deficiency of a Euro per excess cubic meter. In the case of the normal stand with retarded thinning, all excess volumes are small.

**Figure 6.** Capital return rate deficiency per excess volume unit for seven different normal estates, without any restriction in the design of the applied thinning procedure. Only positive values are plotted.

The capital return rate deficiency in terms of Euros per excess volume is shown in Figure 7, as a function of excess volume. There is a restriction of the minimum cutting limit diameter of 238 mm in thinnings from above. Small excess volumes tend toward large values of the observable, since the excess volume is in the denominator. Moderate excess volumes tend to induce a much smaller capital return deficiency than in Figure 6. It is found that increasing the size of trees remaining after thinning is an economical way of increasing expected stand volume, in comparison to extending the rotations (Figures 6 and 7).

**Figure 7.** Capital return rate deficiency per excess volume unit for seven different normal estates, with minimum 238 mm cutting diameter limit in thinnings from above.

The capital return rate deficiency in terms of Euros per excess volume is shown in Figure 8, as a function of excess volume, for the case of heavy proportional thinning. Again, small excess volumes tend toward large values of the observable, since the excess volume is in the denominator. The capital return rate deficiency appears to be greater than in thinnings from above (Figures 6 and 7), for moderate excess volumes somewhat below two Euros per excess volume unit.

**Figure 8.** Capital return rate deficiency per excess volume unit for seven different normal estates, with heavy proportional thinning.

The capital return rate deficiency in terms of Euros per excess volume is shown in Figure 9, as a function of excess volume, for the case of gentle proportional thinning. The capital return rate deficiency appears smaller than in heavy proportional thinning (Figure 8), but greater than in restricted minimum cutting limit diameter (Figure 7). For moderate excess volumes, the gentle proportional thinning in Figure 9 is competitive to merely extending the rotation age in Figure 6.

**Figure 9.** Capital return rate deficiency per excess volume unit for seven different normal estates, with gentle proportional thinning.

#### **4. Discussion**

It appeared from Figure 2 that the optimal rotation age varied significantly. More specifically, the amplitude of variation was 31%. It is of interest to consider whether there are other essential features of the normal stands that would unify. Figure 10 shows the expected value of capital return rate as a function of the state sum of annual volumes per hectare. It is found that at the instant of maximal capital return rate, the state sum of volumes varies significantly, as does the state sum of annual capitalizations (not in the figure). The case with retarded thinning shows more than two times the state sums corresponding to the maximum capital return rate, in comparison to the most fertile stand. The amplitudes of variation of optimal capital return rate and the state sum of annual volumes were 9.8% and 36%, respectively. The amplitude of variation of the optimal state sum of annual capitalizations was 26%.

An obvious reason for the variation in Figure 10 is that sites of different fertility have different time scales in their development. Summing either growing time or some state variables does not unify features unless some time-scaling is applied.

It appears, however, that there is a unifying feature in the seven normal estates. The optimal expected stand volume does not vary much (Figure 11). A consequence is that the excess volumes appearing in Figures 6–9 are additive to an almost unified basic volume, with even the case of retarded thinning differing only slightly. In other words, the amplitude of variation of the optimal expected value of stand volume was 9.1%. This greatly simplifies carbon storage applications.

**Figure 10.** Capital return rate according to Equation (7) for seven different normal estates, without any restriction in the design of the applied thinning procedure, as a function of the state sum of annual volumes.

**Figure 11.** Capital return rate according to Equation (7) for seven different normal estates, without any restriction in the design of the applied thinning procedure, as a function of the expected value of commercial volume per hectare.

We have found from Figures 6–9 that for any of the seven normal estates, the most economical way of increasing carbon storage is to increase the size of trees retained in thinning from above. We even find from Figures 7 and 11 that the economics of carbon storage do not vary widely, even if there are different fertilities and growth rates. Correspondingly, the results of the recent investigation [25] appear to be robust to variations among the sample plots collected in the field.

An iterative search for suitable thinning schedules was done up to the rotation age, giving the greatest capital return rate. After the maximum expected value, any thinning procedure was designed to maximize the expected value of the capital return rate observed within the next five years. Any additional five years was thus considered as a marginal extension of the optimal rotation age.

Within the geographic area, professional practice is not to thin aged spruce stands. This is due to a variety of factors, among them avoidance of wind damage and suspected loss of vigor along with aging. The latter might reduce the recovery of trees remaining after thinning. It is of interest to introduce a boundary condition according to which no thinnings are implemented after the optimal rotation age and then compare it with the results appearing in Figures 2, 3, 6 and 7. It is found from Figures 12 and 13 that excess volumes become greater and the increment of capital return rate deficiency smoother in the absence of post-optimum thinnings.

**Figure 12.** Capital return rate deficiency per excess volume unit for seven different normal estates, without any thinning after the financially optimal rotation age.

**Figure 13.** Capital return rate deficiency per excess volume unit for seven different normal estates, with minimum 238 mm cutting diameter limit in thinnings from above, and no thinning after the financially optimal rotation age.

The capital return rate being strongly affected by the clearcutting price premium, the five-year marginal approach either does or does not trigger further thinnings. If one would adopt another boundary condition, for example considering a 15-year margin, a greater number of thinnings would be triggered. After the optimal rotation age, observables would differ. The expected value of the capital return rate would decline more slowly than in Figures 2–5.

If one would maximize the expected value of the capital return rate observed within any 15-year period after the rotation age corresponding to maximal capital return rates, a chart corresponding to Figure 6 would become meaningless. This is because the average stand volume would decrease with time, and the excess volumes would become negative. Instead, setting a restriction in the cutting diameter limit (as in Figure 7) would yield meaningful results. However, with the cutting diameter limit restriction of 238 mm, the excess volumes are still small and partly negative. With 263 mm cutting diameter limit restriction, the excess volumes become positive. The deficiencies are only slightly greater than in Figure 7 and clearly smaller than in Figure 9. A similar result is not achievable by proportional thinning. If one maximizes the expected value of the capital return rate observed within the next 15 years after the rotation age corresponding to maximal capital return rate, the results become inferior in comparison to Figures 8 and 9.

It is further worth noting that multiple thinnings from above would, along with time, transfer the system from the realm of even-aged forestry towards uneven-aged forestry. Such a transition might reduce vigor [36,37], and such a reduced vigor is not considered in the growth model used in this study.

An applicable carbon rent can be derived from carbon storage market price according to Equation (11). At the time of writing, the market price of carbon dioxide emissions is in the order of 25 Euros per ton. This inserted to Equation (11), together with a 3% discount rate [38–40], results as an annual carbon rent of 0.75 Euros per ton of carbon dioxide.

A cubic meter of commercial trunk volume in the boreal forest stores a ton of carbon dioxide (in living biomass, litter, and soil) [1,2,4,7,41,42]. Correspondingly, one might expect a carbon rent of 0.75 Euros per standing cubic meter. In accordance with a recent study [25], that is enough only up to an excess volume of 30 m3/ha. In other words, the recent result is now confirmed considering variation among sample plots measured in the field. For greater levels of carbon storage, the deficiency of capital return rate per standing cubic meter is greater than the appropriate rent.

Any carbon sequestration trade compensation must be proportional to the carbon storage on the estate level [18]. If it would not be, agents with large carbon inventories would initially suffer heavy and unjustified release taxes. Unlike the carbon trade, there is no particular reason why the carbon rent should be proportional. If the carbon rent is made nonproportional, the marginal rent may reasonably reach or exceed the capital return rate deficiency of Figure 7, also on the right of the figure.

It was recently claimed that changing soil fertility as well as changing temperature sum would change the values of all time derivatives appearing in any growth model, as well as in any carbon storage model [23,25]. Correspondingly, capital return rates, as well as annual capital return deficiencies, would change along with all observables resulting from the integration of time derivatives. This statement is confirmed by Figures 2–5 and 10, where both of the axes contain time-dependent quantities. There is a significant scatter among normal estates corresponding to normal stands of different fertility. On the other hand, quantities not involving time derivatives should be robust to fertility. This is found from the horizontal axes of Figures 6–9 and Figures 11–13, which do not contain time derivatives or their integrals. The vertical axis of the latter figures do contain time-dependent quantities, and correspondingly display significant scatter among the normal estates in relative terms. It is, however, worth noting that the variation of optimal capital return rate in Figures 2, 10 and 11 is surprisingly small. The amplitude was 9.8%, in contrast to 9.1% optimal stand volume. The much bigger variations in other time-dependent quantities may be not only fertility-related but possibly also stand age-related.

In addition to soil fertility, there are other aspects of how the seven normal stands differ. Stem counts and stem size distributions differ, as well as stand age at the time of observation in 2018. Such factors do contribute to the characteristics of the normal estates that the normal stands represent. The results above indicate that for the stand of the oldest age, the first thinning was clearly overdue. Apart from that, the present material does not indicate what would be an optimal stage for the first thinning.

The seven normal stands, representing seven normal estates, were collected from one physical estate in eastern Finland. Correspondingly, the variation among the normal estates reflects variation among stands within the physical estate. The physical estate was completely on mineral soil and relatively productive throughout.

The yield of logs from trees was modeled based on an empirical dataset of 6123 harvested spruce trees. There is no guarantee that such an empirical dataset would apply to all circumstances. The quality of tree trunks varies by stand. Quality requirements of sawlogs vary by sawmill and are subject to change in time. Geographic areas differ, as well as tree species. The yield of logs from trees does differ depending on the type of harvesting (clearcutting, thinning from above/below, continuous cover operations, etc.), stand age, and so on.

Many of the quantitative results in this paper depend on the market prices of roundwood assortments, bare land, silvicultural and harvesting expenses, and carbon emission, as well as the discount rate applied in Equations (10) and (11). Considering eventual changes in the latter two is straightforward in terms of Equations (10) and (11). Changes in the former quantities contribute trivially if they are proportional. In the case of significant nonproportional changes in prices, most of the figures of this paper will have to be redrawn.

#### **5. Conclusions**

The robustness of recent results regarding microeconomics of carbon storage in fertile boreal forest estates is investigated. For all investigated seven normal estates, carbon stocking can be increased at the expense of a capital return rate deficiency. The most economical way of increasing carbon storage is to increase the size of trees retained in thinning from above, in all seven cases.

Normal stands and normal estates show different characteristics in time-dependent features. However, they unify in time-independent characteristics. In particular, the financially optimal expected value of stand volume per hectare does not vary much, which greatly simplifies carbon storage applications.